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PREFACE 


This  report  describes  work  on  new  design  methods  for  filters  used 
in  discrete  data  control  systems.  Design  methods  are  developed  first 
for  sampling  rates  to  minimize  the  bit  requirements  for  each  filter 
coefficient  then  new  design  methods  for  digital  filters  that  minimize 
the  need  for  digital  multiplication  are  described.  Interactive  soft- 
ware for  aiding  design  and  implementation  of  digital  filters  was 
* 

written  and  is  described  in  the  report. 
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INTRODUCTION 


This  report  describes  work  on  new  design  and  implementation 
methods  for  filters  used  in  discrete-data  control  systems. 
Specifically,  the  following  tasks  were  undertaken: 

1.  The  development  of  design  methods  that  use 
the  sampling  interval  as  a design  parameter 
to  minimize  the  bits  required  to  represent 
each  filter  coefficient. 

t 

2.  The  development  of  new  design  methods  for 
digital  control  algorithms  to  minimize  the 
need  for  digital  multiplication. 

I 

3.  The  development  of  interactive  software  to 
aid  in  the  design  and  implementation  of 
digital  control  algorithms. 

4.  A method  of  fault  analysis  for  digital 
control  algorithms. 

Each  task  is  discussed  with  details  given  via  copies  of  each 

report  generated  during  the  contract  period.  These  reports 
are  provided  as  appendices. 

1.  DESIGN  USING  THE  SAMPLING  RATE  AS  A DESIGN  PARAMETER 

Both  digital  and  analog  filter  synthesis  generally  involve 
tradeoffs.  For  instance,  low  ordered  filters  may  have  either 
sharp  rolloff  or  flat  passbands  but  not  both.  High  order  filters 
can  have  excellent  frequency  response  characteristics  but  involve 
a large  number  of  components  or  multiplications,  both  of  which 


increase  errors.  In  sampling  time  synthesis  there  are  tradeoffs 
as  well.  Exact  coefficients  can  be  easily  found  for  quite  a 
few  first  order  filters  but  the  sampling  time  which  yields  such 
coefficients  causes  the  filters  to  have  serious  magnitude  errors 
due  to  aliasing.  On  the  other  hand,  a sampling  time  which  is 
very  short  will  cause  the  filter  frequency  response  to  be  very 
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sensitive  lo  coefficient  quantization.  The  tradeoff  between 
aliasing  errors  and  coefficient  quantization  errors  is  to  be 
kept  in  mind  in  synthesizing  sampling  times  for  bilinear  z- 
transform  filters.  In  fact,  the  tradeoff  considerations  are 
an  important  step  in  the  synthesis  procedure. 

First  Order  Filters 

The  technique  for  synthesizing  sampling  time  is  essentially 
the  same  for  both  first  and  second  order  filters.  However, 
because  there  are  only  two  coefficients  in  the  first  order 
filters  and  because  frequency  independent  bounds  can  be  found 
for  the  first  order  filters,  the  first  order  case  is  developed 
first. 

For  the  design  a realistic  approach  is  to  make  the  magnitude 
(or  phase)  response  errors  as  samll  as  possible  while  retaining 
a short  enough  sampling  interval  to  avoid  aliasing  errors.  One 
means  of  finding  coefficients  which  will  give  small  error  is  to 
generate  a number  of  sets  of  coefficients  and  then  find  the 
truncated  and  bounded  values  for  each.  Next,  take  the  differ- 
ence between  the  designed  coefficients  and  the  respective  quan- 
tized values  and  determine  the  maximum  error  for  each  set  of 
coefficients.  The  set  with  the  smallest  maximum  magnitude  error 
is  then  the  set  of  coefficients  to  use  unless  the  sampling 
interval  associated  with  that  set  is  too  long  to  meet  the  aliasing 
specifications.  A maximum  bound  on  frequency  error  could  be 
given  and  the  sampling  times  and  coefficients  which  give  a mag- 
nitude error  less  than  the  bound  would  be  considered.  If  the 


frequency  response  error  criterion  is  not  met,  then  it  is 
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necessary  to  generate  more  coefficients  using  different  sampling 
times  than  used  previously  and  repeat  the  procedure  above.  If 
the  magnitude  response  error  criterion  is  not  satisfied  after 
several  hundred  sampling  times  have  been  tried  it  would  be 
necessary  to  use  a longer  word  length  to  realize  the  coefficients. 

While  the  procedure  above  seems  quite  _ong,  it  is  possible 
to  combine  all  of  the  steps  of  the  process  into  an  interactive 
computer  program.  The  block  diagram  of  such  a program  is  shown 
in  Figure  1.1. 

The  first  block  of  Figure  1.1  asks  for  input  of  the  analog 

filter  coefficients,  the  word  length  desired,  the  maximum  absolute 

value  for  the  magnitude  response  error,  A | H | , and  the  maximum 

number  of  iterations  to  be  done  before  it  is  decided  that  a 

longer  word  length  is  necessary.  Block  2 initializes  the  sampling 

time  for  a certain  pass.  On  the  first  pass  the  sampling  time 

will  be  set  to  ,01*t  as  an  initial  value  where  t is  the 

max  max 

maximum  value  the  sampling  interval  can  be,  set  to  avoid  aliasing 
errors.  Blocks  3 and  4 are  self-explanatory,  where  block  3 uses 
equations  for  a digital  filter  found  from  an  analog  filter  via 
the  bilinear  z-transform  to  generate  the  digital  coefficients. 

The  fifth  and  sixth  blocks  are  similar  to  each  other.  In  each, 
the  difference  is  found  between  the  designed  (infinite  precision) 
digital  coefficients  and  the  quantized  coefficients.  The  differ- 
ences are  found  for  the  rounded  coefficients  and  then  A | H | is 
derived  by  using  the  magnitude  of  the  desired  and  actual  filters. 
The  same  calculations  are  also  done  for  the  truncated  coefficients 


Then  A ) M | of  the  rounded  values  is  compared  to  the  A | H | of  the 
truncated  values  and  the  smallest  of  those  two  A | H | ’ s is  chosen 
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for  that  sampling  interval.  If  the  A | H | chosen  is  less  than  the 
maximum  allowable  magnitude  error,  specified  in  the  first  block, 
then  the  sampling  time  and  the  associated  digital  coefficients 
are  printed  out  along  with  the  type  of  quantization  to  be  used. 
Also,  if  the  A | H | chosen  on  a particular  iteration  is  smaller 
than  any  chosen  on  any  previous  iteration  then  it  is  stored 
along  with  its  corresponding  sampling  interval  and  the  previously 
stored  values  are  discarded. 

Whether  or  not  A | H | is  less  than  the  previous  smallest 

value,  the  sampling  time  is  increased  by  ,01-t  . If  the  sam- 

max 

pling  time  is  then  less  than  or  equal  to  t a new  set  of  coef- 

max 

ficients,  differences,  and  magnitude  response  errors  is  gener- 
ated. If  the  sampling  time  is  greater  than  t^  then  the 
procedures  of  blocks  15  through  18  are  executed.  If  there  is 
at  least  one  A | H | of  those  considered  which  meets  the  error 
criterion  then  the  sampling  time  which  gave  the  smallest  A | H | 
is  output  along  with  the  A | H | . Otherwise  a test  is  done  to 
see  if  the  maximum  number  of  iterations  have  been  run  through. 

If  they  have  then  it  is  advisable  to  increase  the  word  length 
and  run  through  the  iterations  again.  If  the  maximum  number 
of  iterations  has  not  been  reached  then  a new  set  of  sampling 
times  should  be  tried.  A search  can  be  made  around  the  immediate 
area  of  the  sampling  time  which  gave  the  lowest  a|h|,  using  a 
smaller  sampling  time  increment  for  the  new  iterations.  Another 
possibility  is  to  merely  offset  the  new  sampling  times  from  those 
of  the  previous  pass  by  a certain  amount,  for  example  . 001‘t^^. 

A FORTRAN  program  which  realizes  the  block  diagram  of  Figure  1.1 
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has  been  used  to  illustrate  the  procedure. 

Second  Order  Filters 

The  block  diagram  in  Figure  1.1  for  first  order  filter 
synthesis  servo  as  well  for  second  order  filters.  If  a second 
order  low  pass  section  is  being  designed,  a possible  evaluation 
technique  would  be  to  evaluate  A | H | at  SI  = 0 and  at  the  3-db 
point  for  every  sampling  time  considered.  For  a bandpass  struc- 
ture A | H | could  be  calculated  at  the  3-db  points  and  the  pole 
frequency.  In  the  example  programs  the  user  is  allowed  to 
choose  what  radian  frequencies  the  magnitude  error  is  to  be 
evaluated  at,  and  how  many  frequencies  are  to  be  evaluated. 

The  block  diagram  in  Figure  1.1  allows  every  sampling 
interval  and  corresponding  set  of  coefficients  which  have  a 
A | H | smaller  than  the  maximum  error  bound  to  be  printed  out. 

The  reason  for  this  is  very  simple.  In  some  cases  a certain 
sampling  rate  will  be  more  desirable  than  another  even  if  it 
does  not  give  the  minimum  magnitude  response  error.  Such  a 
case  occurs  when  the  clock  rate  for  the  filter  is  limited  to  a 
certain  range  of  values.  By  printing  all  sampling  times  which 
have  magnitude  errors  within  the  desired  bound,  there  is  more 
design  freedom  permitted. 

So  far  the  discussion  has  centered  about  magnitude  response 
design.  However,  a bilinear  z-transform  digital  filter  will 
not  generally  have  the  same  phase  response  as  the  corresponding 
analog  filter.  However,  digital  filters  do  have  phase  response 
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and  it  is  sometimes  desirable  to  retain  the  accuracy  of  that 
response.  It  is  possible  to  make  the  phase  response  accurate 
in  the  same  manner  that  was  used  for  the  magnitude  response. 

In  fact,  the  block  diagram  of  Figure  1.1  can  be  used  for  the 
phase  by  replacing  |h|  by  0 in  the  diagram.  A longer  flow 
chart  could  be  developed  which  allows  filters  with  accurate 
phase  and  magnitude  to  be  designed. 

In  terms  of  design  limitations,  the  primary  considerations 
are  those  of  aliasing  and  processor  (or  component)  speed.  When 
deriving  a program  to  find  a sampling  time  which  results  in 
small  frequency  response  error,  it  is  necessary  to  put  an  upper 
bound  on  the  length  of  the  sampling  interval  to  reduce  aliasing 
errors.  If  the  sampling  rate  is  too  low,  a filter  will  be 
aliased  to  the  point  where  it  no  longer  performs  its  designed 
task.  To  reduce  aliasing  errors,  it  could  be  required  that  the 
sampling  rate  be  at  least  ten  times  the  highest  pole  frequency. 
There  are  various  such  rules  of  thumb  aimed  at  avoiding  aliasing 
errors  and  which  ever  is  appropriate  should  be  used. 

Digital  hardware  is  limited  in  terms  of  clock  rates  it  can 
operate  at.  Some  computers  can  perform  an  instruction  in  a 
matter  of  nano-seconds  while  others  require  several  micro-seconds 
to  do  the  same  instruction.  Similarly,  discrete  digital  com- 
ponents such  as  multipliers,  adders,  and  shift  registers  are 
limited  in  speed.  When  designing  a digital  filter  it  is  important 
to  realize  the  constraint  of  digital  hardware  speed  on  sampling 
time.  If  a computer  program  is  used  to  realize  a filter,  the 


program  may  be  ten,  twenty,  or  even  over  a hundred  instructions 


long.  Often,  large  filter  structures  are  better  realized  using 
discrete  digital  components  which  are  dedicated  to  the  filtering 
because  the  discrete  components  have  an  advantage  in  speed  over 
a computer  program.  Whatever  method  of  realization  is  used 
though,  tt  2 design  should  not  allow  the  clock  rate  of  the  filter 
exceed  the  speed  of  the  structure  used  to  realize  it. 

Example  Designs 

A commonly  used  filter  design  is  the  maximally  flat,  or 
Butterworth,  filter.  The  low  pass  Butterworth  filter  has  the 
property  that  the  filter  magnitude  response  is  as  flat  as  possible 
at  co  = 0.  For  the  present  example,  a fifth  order  Butterworth 
low  pass  digital  filter  is  synthesized  using  the  method  de- 
scribed in  the  first  part  of  this  chapter.  The  analog  transfer 
function,  H(s),  is  given  by 

h(s)  * s+r  ~~i 2 1 a-u 

s +.6l8s+1.0000  s +1.618s+1.0000 
so  that  H(s)  has  unity  gain  and  unity  bandwidth.  The  example 
demonstrates  the  use  of  both  the  first  and  second  order  synthesis 
programs.  The  bilinear  z-transform  allows  the  transfer  function 
to  be  broken  up  into  first  and  second  order  cascade  sections  so 
there  is  no  partial  fraction  expansion  to  worry  about.  For  the 
example,  the  assumptions  are  that  an  eight  bit  word  and  a magnitude 
error  of  less  than  10  ^ are  desired. 

The  first  order  section,  H^(s),  of  H(s)  is  given  by 
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Figure  1.2  shows  the  sequence  of  interactive  inputs  to  the  program. 
The  first  three  inputs  are  self  explanatory.  After  the  word 
length  was  input,  the  program  used  100  sampling  times  between  0 
and  the  maximum  sampling  time  allowed.  An  appropriate  sampling 
time  was  not  found  and  the  graph  of  Figure  1.3  resulted.  The 
graphs  are  not  meant  to  be  an  absolute  means  of  measuring  the 
error  of  the  filter  but  merely  a way  of  determining  whether  to 
proceed  or  to  try  another  word  length  or  error  bound . The  next 
input  given  in  Figure  1.2  was  a l(one)  to  indicate  that  on  the 
next  iteration  the  same  range  of  sampling  times  was  to  be  used 
but  the  sampling  times  would  be  offset  from  the  previous  set  of 
sampling  times.  The  amount  that  the  second  set  was  offset  was 
one-tenth  of  the  spacing  of  the  first  set  of  sampling  times. 
Therefore,  there  was  a sampling  time  selected  between  each  of 
the  first  sampling  times.  Again  the  error  bound  was  not  met, 
resulting  in  Figure  1.4  which  is  very  similar  to  Figure  1.3. 

Rather  than  continue  on  the  same  track,  it  was  felt  that  a better 
approach  would  be  to  "blow  up"  the  region  around  the  sampling 
interval  which  gave  the  minimum  error.  100  sampling  times  were 
chosen  between  the  two  sampling  intervals  which  were  adjacent  to 
the  point  which  gave  the  least  error.  The  input  of  2 in  Figure 
1.2  resulted  in  the  expansion  about  the  minimum  point  and  the 
graph  of  Figure  1.5.  The  error  bound  was  still  not  met  but  there 
seemed  to  be  promise  so  another  expansion  was  done.  Figure  1.7 
shows  that  the  error  bound  was  finally  satisfied  by  three  sampling 
times.  The  output  lists  the  three  sampling  times  which  allowed 
the  error  criterion  to  be  met  and  the  corresponding  coefficients 


and  the  type  of  coefficient  quantization  to  be  used  on  each  set 
of  coefficients.  Figure  1.6  shows  the  graph  of  the  final  expansion. 
The  listing  on  Figure  1.7  prints,  as  a final  set  of  values,  the 
minimum  magnitude  error  found  and  the  sampling  time  which  gave 
the  minimum  error.  If  there  had  been  no  sampling  times  which 
caused  the  error  bound  to  be  satisfied  on  that  last  round,  then 
it  probably  would  have  been  necessary  to  use  a longer  word  length 
or  accept  a slightly  relaxed  error  bound.  Sometimes  it  is  possible 
to  do  enough  passes  to  find  a sampling  time  which  gives  low 
enough  error  but  in  order  to  set  the  sampling  time  it  would  require 
an  infinitesimal  adjustment  and  so  the  sampling  time  is  not  prac- 
tically realizable.  Even  with  programmable  clocks,  the  adjustment 
is  usually  only  down  to  about  10  ^ so  adjustments  below  that 
level  are  not  possible. 

The  two  second  order  sections,  H^Cs)  and  H^(s),  are  given 
by 

H (S)  = (1.3) 

s +.618s+1.0000 

Ho(s)  = — p (1.4) 

s +1.6l8s+1.0000 

Since  there  is  no  one  general  frequency  which  results  in  a maximum 
for  the  partial  derivatives  in  the  expansion  AH  [1,2],  then 
several  frequencies  should  be  chosen  to  check  those  partial 
derivatives.  In  the  examples,  four  frequencies  were  chosen  for 
each  section.  Three  were  in  the  past.oand  and  one  was  in  the 
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transition  band  of  each  filter  section.  The  graphs  of  Figures 
1.8  and  1.9  show  the  relative  errors  of  the  filters  against  the 
sampling  times  for  H^Cz)  and  H^(z),  respectively.  Some  of  the 
errors  are  so  large  that  most  of  the  error  points  appear  to  be 
very  small  but,  in  reality,  only  a few  of  the  plotted  points 
resulted  in  filters  which  satisfied  the  error  bounds. 

The  second  order  filter  errors  behave  much  differently  than 
the  first  order  filter  errors.  The  first  order  errors  tend  to 
decrease  as  the  sampling  interval  gets  longer,  while  the  second 
order  errors  tend  to  increase.  Also,  the  second  order  errors, 
with  a few  exceptions  as  noted  on  the  graphs,  are  generally 
much  lower  than  the  first  order  errors.  Therefore,  it  seems 
that  the  word  length  restrictions  on  a filter  would  come  from 
the  filter's  first  order  section  or  sections.  The  design  tech- 
nique, then,  should  state  that  the  first  order  sections  should 
be  synthesized  before  the  second  order  sections  in  order  to  get 
a good  bound  on  the  word  length  requirements. 

After  considerable  effort  with  more  examples  this  approach 
appeared  to  be  somewhat  limited,  in  fact,  the  method  does  not 
generally  work.  To  simplify  the  coefficient  problem  and  to 
totally  eliminate  the  multiplier,  a different  procedure  was  tried 
as  shown  in  the  next  section. 
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2.  NEW  DESIGN  METHODS  FOR  DIGITAL  CONTROL  ALGORITHMS 

The  use  of  conventional  digital  control  algorithms  requires 
the  implementation  of  a recursive  difference  equation  on  a com- 
puter. Here  multiplication  and  addition  operations  are  employed. 
The  finite  length  arithmetic  causes  roundoff  errors  to  occur  and 
these  errors  can  be  highly  dependent  on  the  sampling  interval. 

Of  particular  importance  are  the  coefficient  rounding  errors, 
because  if  the  coefficients  are  not  rounded  properly  the  al- 
gorithm may  be  unstable  or  not  exhibit  the  desired  magnitude 
and  phase  characteristics.  Filter  designs  employing  the  bilinear 
z-transform  are  very  subject  to  these  errors  because  the  coef- 
ficient sensitivity  becomes  increasingly  critical  as  the  filter 
order  increases. 

To  avoid  the  whose  process  of  multiplication  and  to  decrease 
the  filter  sensitivity  to  the  sampling  interval,  a radical  new 
approach  to  filter  implementation  has  been  investigated.  The 
initial  idea  is  that  any  filter  can  be  approximated  by  a non- 
recursive filter  with  a transfer  function  of  the  form 

N -k 

H(z)  = E a z 

k=0 

Now  if  the  filter  can  be  structured  in  such  a way  to  make  the  coef- 
ficients powers  of  2,  all  multiplication  can  be  done  via  shifting. 
If  this  is  done,  high  order  nonrecursive  structures  can  be  used 
to  approximate  even  low  order  recursive  filters  with  definite 
advantages. 

The  first  advantage  is  that  implementation  can  be  done 
directly  using  large  scale  integrated  circuits.  Such  an  imple- 
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mentation  allows  the  shifting  operations  to  be  designated  under 
computer  control  so  that  a design  can  be  used  to  implement  a 
variety  of  different  filters. 

The  second  advantage  is  that  linear  phase  can  be  implemented 
and  this  can  be  very  useful  in  control  system  design. 

The  details  of  this  work  are  provided  by  the  report  in 
Appendix  1.  Here  the  design  methodology  along  with  several 
examples  provide  the  necessary  background  for  this  implementation. 

Experiments  have  also  been  tried  on  recursive  filters  using 
the  implementation  procedure  given  in  the  report.  The  results 
up  to  this  time  are  mixed,  however,  using  the  integer  design 
approach  we  feel  that  good  recursive  filter  implementations  can 
be  obtained. 
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3.  INTERACTIVE  SOFTWARE 

Four  interactive  digital  filter  design  packages  were  written 
that  are  valuable  to  the  designer  of  digital  control  algorithms 
and/or  digital  signal  processing  algorithms.  These  are: 

a.  A Digital  Filter  Design  Program  Utilizing 
the  Bilinear  Z Transform 

b.  Programs  for  Weighted  Least  Squares  Design 
of  Nonrecursive  and  Recursive  Digital 
Filters 

c.  A Fortran  IV  Design  Program  for  Low-Pass 
Butterworth  and  Chebychev  Digital  Filters 

d.  A Fortran  IV  Design  Program  for  Butterworth 
and  Chebychev  Band-Pass  and  Band-Stop 
Digital  Filters 

All  of  the  programs  are  written  in  FORTRAN  and  run  on  the 
DEC  PDP-11  and  the  DEC  GT-40  Graphics  System.  The  detailed 
descriptions  are  given  in  Appendices  B,  C,  D and  F.  Card 
decks  were  sent  to  WPAFB  during  the  course  of  the  contract. 

4.  FAULT  DETECTION 

Because  of  the  need  to  understand  how  well  a digital  control 
algorithm  is  operating,  some  work  was  done  on  detecting  when  a 
discrete  time  algorithm  is  not  operating  correctly  due  to  hard- 
ware or  software  failure.  A parameter  identification  algorithm 
was  used  with  the  method  being  described  via  the  report  of 
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ABSTRACT 


The  difference  routing  digital  filter  (DRDF)  is  a FIR 
filter  whose  coefficients  are  equal  to  zero,  or  integral  powers 
of  two.  The  basic  DRDF  structure  is  reviewed,  and  two  coefficient 
restrictions  are  detailed  that  will  insure  bounded  input,  bounded 
output  stability  as  well  as  a finite  impulse  response.  Next, 
three  parallel  structures  are  presented.  Each  of  these  new 
structures  will  significantly  reduce  the  RMS  error  between  the 
desired  impulse  response  and  the  actual  filter  response.  The 
optimum  structure  appears  to  be  a filter  with  a parallel  struc- 
tured transversal  part  with  integer  valued  taps  followed  by  a 
recursive  part  that  in  the  low  pass  case  is  a digital  integrator. 
For  this  new  structure,  an  analysis  is  given  of  the  RMS  error 
performance  in  both  the  time  and  frequency  domain.  This  analysis 
is  supported  by  extensive  computer  simulation  results. 
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INTRODUCTION 


Finite  impulse  response  (FIR)  digital  filter  structures  are 
attractive  in  a variety  of  applications.  Among  their  advantages 
are  the  inherent  stability  and  the  ease  of  realizing  a linear 
phase  characteristic.  Numerous  methods  now  exist  [ 1 ] — [ 3 ] for 
the  design  of  FIR  structures. 

One  disadvantage  of  conventional  digital  FIR  filters,  in 
many  applications,  is  the  slow  operating  speed  due  to  the  large 
number  of  required  multiplies.  Various  methods  [4],  [5]  and 
[6]  have  been  proposed  in  the  past  to  reduce  or  eliminate  this 
multiplier  requirement.  This  paper  focuses  on  the  low  pass 
difference  routing  digital  filter  (DRDF)  [4].  This  filter 
structure  consists  of  a transversal  part  with  coefficients 
restricted  to  be  zero,  or  integral  powers  of  two.  As  originally 
proposed,  the  DRDF  is  limited  in  the  minimum  RMS  error  that  can 
be  obtained  between  the  ideal  and  actual  filter  impulse  response. 

To  reduce  this  error,  three  parallel  structures  are  presented, 
each  of  which  can  significantly  reduce  the  RMS  error  between  the 
desired  filter  response  and  the  actual  filter  response.  These 
three  methods  are  all  structurely  similar,  but  there  are  distinct 
differences  in  the  design  philosophy  used.  In  the  first  two 
methods,  a parallel  structure  is  created  that  approximates  the 
error  that  would  have  occurred  in  the  original  design.  This 
error  signal  is  added  in  such  a way  as  to  provide  overall  error 
reduction.  In  both  cases,  the  parallel  structure  can  be  imple- 
mented with  minimal  additional  hardware.  The  third  approach, 
which  appears  to  be  optimum,  also  uses  a parallel  filter  structure. 
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In  this  latter  case,  the  purpose  of  the  parallel  structure  is 
to  provide  integer  valued  taps  that  more  accurately  approximate 
the  desired  filter  function.  This  third  approach  has  two  prin- 
ciple advantages  over  the  first  two  methods:  (1)  reduced  hard- 
ware requirements,  and  (2)  a more  straightforward  and  logical 
design  procedure. 

Several  examples  of  improved  performance  with  the  new 
structures  are  given  using  computer  simulation.  An  analysis 
of  the  RMS  error  performance  of  the  optimum  structure  is  also 
given  in  both  the  time  and  frequency  domain.  The  analysis  is 
seen  to  agree  very  favorably  with  computer  simulation  results. 

This  work  was  all  done  with  low  pass  filter  designs. 
However,  similar  results  can  also  be  obtained  for  both  band 
pass  and  high  pass  circuits. 


j 
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II.  LOW  PASS  DRDF  STRUCTURE 


The  structure  of  a low  pass  DRDF  is  shown  in  Figure  1.  The 

coefficients  a ....  a.,  , of  the  transversal  part  are  restricted 
o N-l 

to  be  zero  or  integer  powers  of  two.  The  recursive  part  is  a 
digital  integrator  with  a single  coefficient  equal  to  minus 
one.  The  coefficients  of  the  transversal  part  are  chosen  to  be 
approximations  to  the  differences  between  successive  values  of 
the  desired  filter  impulse  response  li^[nT] 

a * hD[jT]  - hD[ (j-l)T]  (1) 

Thus,  if  the  input  signal  to  the  filter  is  a unit  impulse. 


6[nT]  = 1 
= 0 


for  n = 0 
for  n ^ 0. 


(2) 


Therefore,  the  output  y(nT)  will  be  an  approximation  to  the 
desired  impulse  response  itself 

N-l 

h_[nT]  * y[nT]  = y[  (n-l)Tl  + Z a.  6[(n-j)T]  (3) 

U j=0  J 

The  use  of  a digital  integrator  places  one  restriction 
on  the  a^  coefficients  to  insure  a finite  impulse  response 
(FIR)  filter,  the  sum  of  the  a^  coefficients  must  be  zero. 

This  is  shown  as  follows.  If  the  filter  is  FIR,  then  we 
desire  y[(N-l)T]  to  be  zero. 

Hence : 

y [ (N-l)T  ] = y [ (N~2)T ] + a^  = 0 (4) 

but,  N-2 

y[(N-2)T]  = Z a. 

j=0  J 


(5) 
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Therefore: 


N-l 


E a. 

J-o  J 


0 


(6) 


The  restriction  of  equation  six  is  not  a practical  problem. 
Since  the  desired  finite  impulse  response  in  practice  will 
always  damp  out  to  zero  at  (N-l)T  as  in  Figure  2,  we  have. 

N-l 

h [ (N-l)T]  = 0 = y[(N-l)T]  = E a (7) 

u j=0  J 


Therefore,  the  restrictions  on  the  a_.  coefficients  are: 
a = 0,  +2K  for  K = 0,  1,  2,  ... 

j 

N-l 

E a.  = 0 


(8) 
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III.  DRDF  OPERATION  AND  DESIGN 

The  design  of  a DRDF  is  based  on  approximating  a desired 
finite  impulse  response.  Consider  a desired  impulse  response 
h [nT],  the  first  few  samples  of  which  are  shown  in  Figure  3. 
Without  loss  of  generality,  assume  hQ[0]  is  zero.  Further, 
consider  that  h^nT]  has  been  amplitude  scaled  (h^lnT]  = F-h^JnT]) 


such  that  the  maximum  change  is: 


Era 


max  I hg [ jT J - h^ [ ( j — 1 > T ] | = 2 


where  Km  is  the  largest  exponent  being  considered  in  the  design. 
Thus,  the  DRDF  will  approximate  a scaled  version  of  the  desired 
impulse.  This  scaled  value  is  then  multiplied  by  a scale  factor 
(F)  to  give  the  desired  impulse  response  as  in  Figure  4.  This 
will  insure  that  the  coefficient  values  will  be  able  to  follow 
the  maximum  slope  of  the  scaled  impulse  response,  hg[nT]. 


Note  that: 


Km 


max  |hD[jT]  - hD[(j-l)T] | = F 2 


Therefore,  as  Km  increases,  the  scale  factor  F will  get  smaller. 

The  value  of  the  first  coefficient,  a is  selected  from  0,  +2  ; 

o - 

K = 0,  1,  . . . , Km  so  as  to  be  closest  to  the  first  change 
hg[T)  - hg[0].  Since  hg[0]  equals  zero,  we  have 


a - h [T]  (9) 

o s 


The  design  proceeds  recursively,  selecting  a^  from  0,  +1,  +2 
to  minimize: 

j-1 

| a.  - {h  [ (j+l)T ] -l  a } | 

J S k=0  K 


(10) 
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Figure  5 is  a comparison  of  the  entire  desired  impulse  response 
sequence  hQ[nT]  and  the  error  for  the  DRDF  approximation  for 
Km=4.  Figure  6 is  a comparison  of  the  desired  magnitude  response 
and  the  error  for  the  DRDF  magnitude  response.  For  this  example 
T = 0.05  sec  and  N = 200.  The  "IDEAL"  Chebychev  impulse  response 
used  in  this  and  all  subsequent  examples  was  obtained  from  syn- 
thetic division  of  the  H(Z)  found  using  the  impulse  invariant 
design  [7]. 

Figure  7 is  a plot  of  the  RMS  error  between  the  desired 
impulse  response  and  the  actual  DRDF  impulse  response  for  the 
four  pole  Chebychev  filter.  The  RMS  error  is  expressed  as  a 
percentage  of  the  peak  value  of  the  impulse  response,  and  it 
is  plotted  vs  Km.  It  is  seen  that  little  improvement  results 
beyond  a Km  of  3 or  4.  The  parallel  structure  introduced  in  the 
next  section  provides  a method  of  significantly  reducing  this 
RMS  error. 

IV.  PARALLEL  DRDF  STRUCTURE 

It  was  shown  in  the  previous  section  that  increasing  Km 
beyond  3 or  4 dcc?=  little  to  further  reduce  the  percentage  RMS 
error.  There  may  be  many  applications  where  further  improve- 
ment is  desirable.  One  way  to  do  this  is  to  both  double  the 
sampling  rate  and  also  the  number  of  transversal  stages.  For 
example,  doubling  the  sampling  rate  and  doubling  the  number  of 

transversal  stages  will  cut  the  RMS  error  in  half.  Since  it 
may  not  always  be  possible  to  double  the  sampling  rate,  and 

since  doubling  the  required  number  of  stages  is  not  attractive, 
another  alternative  is  desirable.  Three  different  alternative 


designs  are  considered  below. 
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In  the  first  method,  an  error  sequence  e[nT]  is  defined  as 
the  difference  between  the  desired  impulse  response  and  that 
actually  generated  by  the  DRDF.  That  is: 

e [nT  1 = hulnT]  - hJnT]  (11) 

If  the  error  sequence  e[nT]  of  equation  11  could  somehow 
be  approximated,  e[nT]  and  added  to  the  DRDF  output,  the  new 
signal  h^[nT]  = h^[nT]  + e[nT]  would  be  a better  approximation 
to  the  desired  signal.  Since  e[nT]  is  itself  a finite  duration 
sequence,  it  is  possible  to  approximate  it  with  a second  DRDF 
filter  as  shown  in  Figure  8.  These  two  parallel  filters  can 
share  much  of  the  basic  DRDF  hardware  as  shown  in  Figure  9. 

Note  that  the  parallel  DRDF  will  have  its  own  scale  factor  F^. 

In  some  cases,  F^  can  itself  be  satisfactorily  approximated 
by  an  integral  power  of  2,  however,  in  general  this  is  not  the 
case . 

Conceptionally , this  process  of  approximating  the  DRDF 
error  could  be  continued  to  two,  three  or  more  parallel  stages. 
Of  course,  at  some  point  it  will  be  more  expedient  to  use  a 
conventional  filter  structure. 

Figure  10  is  a plot  of  the  percentage  RMS  error  vs  Km  for 
the  basic  DRDF  and  for  one  and  two  parallel  stages.  This  is  for 
the  Chebychev  filter  used  in  previous  examples.  It  is  seen  that 
the  RMS  error  is  reduced  by  a factor  of  about  3 each  time  a 
parallel  branch  is  added.  Thus,  RMS  errors  well  below  1%  of 
the  peak  value  of  the  impulse  response  are  feasible  with  this 


approach. 
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The  parallel  filter  does  not  approximate  the  error  wave- 
form nearly  as  well  as  the  basic  DRDF  matches  the  original 
desired  impulse  response.  This  is  because  the  error  sequence 
is  quite  noise-like  with  rapid  changes.  The  error  waveform  for 
the  Chebychev  filter  was  shown  in  Figure  5.  Except  for  the 
final  sequence  values,  the  error  signal  is  very  much  like  white 
noise.  The  sinusoidal  appearance  of  the  final  sequence  values 
is  due  to  the  fact  that  the  small  ripple  values  in  the  desired 
impulse  response  is  being  matched  by  a zero  output  from  the 
DRDF. 

It  has  been  found  that  a consistently  better  approximation 
to  the  error  signal  may  be  made  by  roughly  quantizing  the 
error  signal  to  integer  powers  of  2 or  zero.  Thus,  the  new 
filter  shown  in  Figure  11  would  be  similar  to  that  of  Figure  9, 
but  without  a second  integrator.  This  is  the  second  design 
method . 

Figure  12  is  a plot  of  the  percentage  RMS  error  vs  Km  for 
the  basic  DRDF  and  for  one  and  two  parallel  stages  where  the 
parallel  sections  are  rough  quantizations  of  the  error  signals 
to  integer  powers  of  2 or  zero.  A comparison  with  Figure  10 
shows  that  this  second  approach  is  clearly  better.  Similar 
improvements  for  other  filters  have  also  been  found,  and  the 
results  of  Figures  10  and  12  may  be  considered  typical. 

The  third  method  is  aimed  directly  at  the  reason  why  the 
basic  DRDF  error  performance  does  not  improve  as  Km  goes  beyond 
3 or  4.  This  is  because  as  Km  increases,  the  allowable  tap 
values  are  spread  further  apart.  If,  however,  as  Km  increased, 
all  the  integer  values  were  allowable,  then  clearly  the  quanti- 
zation would  improve  and  result  in  reduced  RMS  error.  It  is 


worth  noting  at  this  point  that  uniform  quantization  intervals 
of  any  desired  value  could  be  achieved  with  appropriate  scaling. 
This  then  represents  a rough  quantization  of  the  transversal 
coefficients  with  the  subsequent  integration  acting  to  smooth 
the  overall  impulse  response. 

Because  of  the  speed  and  cost  benefits,  it  is  desirable  to 
obtain  integer  value  taps  with  shifting  and  adding  rather  than 
by  the  use  of  hardware  multipliers.  A direct  approach  would 
be  to  have  a parallel  filter  section  for  various  integer  values 
of  2.  In  this  case,  each  tap  weight  would  be  constructed  from 
its  binary  equivalent.  For  example,  for  a tap  value  of  5 (101) 
there  would  be  three  parallel  filter  sections  with  connections 
made  to  the  first  and  third,  but  not  to  the  second  section. 

Each  section  would  have  its  own  adder.  Figure  13  shows  such 
a filter  which  is  capable  of  producing  tap  values  0 to  7,  and 
with  suitable  two's  complement  circuitry,  -7  to  +7. 

A similar  approach,  but  one  with  further  hardware  economies, 
is  to  permit  both  positive  and  negative  values  of  the  integer 
values  of  two.  Refer  to  Table  I.  This  shows  how  the  integers 
up  to  42  could  be  implemented  with  just  three  parallel  sections. 
The  first  section  could  have  up  to  five  shifts,  the  second  sec- 
tion up  to  three  shifts,  and  the  third  section  a single  shift. 
Thus,  31  is  implemented  as  32-1  rather  than  as  16+8+4+2+1. 

The  optimum  implementation  of  this  concept  will  be  dependent 
on  the  application  and  device  technology.  One  possible  struc- 
ture is  presented  by  Kishi  et  al.  [8]. 

Thus,  we  have  considered  a third  method  of  reducing  RMS 
error  that  of  creating  integer  value  taps.  We  have  also  looked 
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at  several  possible  implementations  of  the  integer  tap  concept. 
A major  advantage  of  the  integer  tap  approach  is  that  separate 
scale  factors  (hardware  multipliers)  are  not  required  in  the 
integer  tap  approach.  A second  advantage  with  this  method  is 
that  error  performance  continues  to  improve  as  the  number  of 
shifts  is  increased.  In  the  basic  DRDF  and  the  other  parallel 
filter  approaches,  improvement  leveled  off  beyond  a Km  of  3 or 
4. 

Figure  14  compares  the  performance  of  the  DRDF  with  the 
integer  tap  approach  for  the  four  pole  Chebychev  filter.  In 
the  graph  at  Km  = 4,  the  basic  DRDF  is  allowing  taps  values 
of  0,  +1,  +2,  +4,  +8  and  +16,  while  the  integer  approach  is 
allowing  all  the  integer  values  0,  +1,  +2,  +3  . . . . +16. 

For  this  filter,  the  integer  approach  surpasses  the  best  re- 
sults of  the  other  methods  at  a maximum  integer  value  of  +138. 

Figure  14  compares  the  performance  of  the  basic  DRDF  with 
the  integer  tap  approach  for  the  four  pole  Chebychev  filter. 

The  design  of  the  integer  taps  would  proceed  in  a recur- 
sive fashion  similar  to  the  basic  DRDF  design.  Tap  values  a^ 
would  be  selected  from  the  allowable  integers  0,  +1,  +2,  ... 
+MAX  to  minimize: 


j-1 

- {h  t(j+l)T]  - Z a } | 
k=0  K 


a . 
J 


(12) 


V.  TIME  DOMAIN  ERROR  PERFORMANCE 


An  important  measure  of  performance  will  be  the  RMS  error 
between  the  desired  impulse  response  sequence  h^[nT]  and  the 
actual  sequence  h^tnT].  In  the  time  domain,  the  RMS  error  can 
easily  be  calculated  for  any  filter  design  from: 


' N-l  h [nT]  - hA[nT) 

l — 


n=0 


N 


(13) 


where  the  subscript  T stands  for  time-domain. 

It  is  desirable  to  be  able  to  estimate  what  this  error  may 

be  for  a particular  filter  without  going  through  the  actual 

design  procedure.  In  this  section,  the  time-domain  RMS  error 

of  the  integer  tap  approach  is  estimated. 

The  errors  measured  at  the  filter  output  may  be  assumed  to 

2 

be  uniformly  distributed  with  zero  mean  and  variance  Q / 1 2 [8]. 

In  the  case  of  integer  taps  Q = 1.  Since  the  mean  is  zero,  the 
RMS  error  will  be  l//l2. 

In  actual  practice,  the  designer  would  want  to  know  how 
the  RMS  error  compared  with  the  peak  value  of  the  impulse  response 
sequence.  Therefore,  it  is  desirable  to  have  an  estimate  of  the 
peak  value  of  the  impulse  response  sequence.  This  estimate  can 
be  made  using  the  following  rules. 

1.  The  main  lobe  of  a typical  high  order  low  pass  filter  will 
have  a width  that  is  approximately  equal  to  the  reciprocal  of 
the  cutoff  frequency. 

2.  The  average  slope  of  the  main  lobe  will  be  about  one  half 


its  maximum  value. 
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Th'-refore,  an  estimate  of  the  peak  output  of  a DRDF  can 
be  made  given  the  filter  cutoff  frequency,  fco,  the  sampling 
interval,  T,  and  the  maximum  integer,  Im. 


= Im  . 1 

2T  ’ 2 fco 


(14) 


In  Table  II,  a comparison  is  made  of  actual  and  estimated 
RMS  errors  as  a percentage  of  the  peak  value  of  the  impulse 
response  sequence.  The  percentage  estimates  are  obtained  by 
dividing  the  RMS  estimates  obtained  from  equation  (13)  by  the 
peak  value  estimates  from  equation  (14).  In  all  cases,  the 
estimates  represent  a conservative  bound  on  the  actual  error. 


The  RMS  error  estimate  as  a percentage  of  the  peak  value  (E') 
is  thus  seen  to  be: 


(15) 


Therefore,  for  a fixed  sample  rate,  the  percentage  RMS  error  is 


inversely  proportional  to  the  maximum  integer  value,  Im. 


VI.  FREQUENCY  DOMAIN  ERROR  PERFORMANCE 

Since  filter  performance  requirements  are  often  given  in 
terms  of  the  frequency  domain,  it  is  important  to  evaluate  the 
frequency  response  error  performance. 

The  error  sequence  e[nT]  is  defined  in  equation  (11)  as 
the  difference  between  the  desired  impulse  response  and  that 
actually  generated  by  the  filter.  In  the  frequency  domain,  the 
magnitude  of  the  error  at  any  frequency  ui^  can  be  obtained  by 
evaluating  the  z transform  of  e[nT]  at  z = eJWiT.  Therefore, 


the  magnitude  response  of  the  error  may  be  written  as: 


A t T N_1 

E(u>)  S E(z=eJa")  = E , -nT 

q e l nT  J z (16) 

If  equation  (15)  is  evaluated  for  a set  of  frequencies  wo> 
u^,  ...  u^_^,  then  we  can  define  an  expression  for  t he  RMS  error 
in  the  frequency  domain. 


We  can  see  from  equation  (15)  that  the  value  of  E(u)  and 
hence,  the  value  of  the  RMS  error  Eu  is  a function  of  the  number 
of  coefficients,  N in  the  FIR  filter  used  to  generate  the  impulse 
response.  In  1973,  Chan  and  Rabiner  [9]  showed  that  Eu,  the  RMS 
error  in  the  frequency  domain  is  found  from: 

Eu  = /N  Et  (18) 

here  ET  is  the  time  domain  RMS  error  defined  in  equation  (13). 

Estimates  of  the  frequency  domain  RMS  error  can  be  derived 
from  the  time  domain  RMS  error  estimate  by  application  of 
equation  (18).  Table  III  is  a comparison  of  the  estimated 
frequency  response  errors  with  the  actual  errors  for  the  Cheby- 
chev  filter  previously  used. 

VII.  OPERATION  EXAMPLES 

Specific  examples  are  given  in  this  section  of  the  time  and 
frequency  domain  performance  of  DRDF  structures  compared  with 
the  ideal  time  and  frequency  responses.  The  marked  improvement 
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using  the  integer  filter  is  also  shown.  Because  the  performance 
is  so  close  to  the  ideal,  especially  if  an  integer  filter  is  used, 
only  the  error  (difference  from  ideal)  is  plotted  along  with  the 
ideal  waveforms. 

Figure  15  shows  the  ideal  impulse  response  for  the  Cheby- 
chev  filter  used  in  previous  examples.  Also,  plotted  with  an 
expanded  amplitude  scale  are  the  error  waveforms  for  an  integer 
filter  with  Im  = 16.  The  improvement  gained  with  the  integer 
filter  is  readily  apparent  by  comparing  Figure  15  with  the  basic 
DRDF  results  of  Figure  5.  There  are  reductions  in  both  the  peak 
and  RMS  errors.  This  same  improvement  is  mirrored  in  the  fre- 
quency domain  as  seen  in  Figure  16  which  shows  the  ideal  magni- 
tude response  along  with  the  magnitude  errors  for  the  integer 
filter.  Compare  Figure  16  with  the  results  shown  in  Figure  6. 

Figure  17  shows  the  ideal  step  response  for  the  same  Cheby- 
chev  filter.  Again,  the  error  waveforms  are  plotted  on  the  same 
time  scale  and  we  see  the  improvement  achieved  with  the  use  of 
the  integer  filter.  The  step  response  also  indicates  that  the 
filter  is  BIBO  stable.  Very  similar  results  are  obtained  with 
all  other  low  pass  filter  designs. 

VIII.  SUMMARY 

The  structure  and  performance  of  the  difference  routing 
digital  filter  (DRDF)  has  been  explored.  Design  restrictions 
and  the  basic  filter  design  have  been  detailed.  Examples  of 
the  RMS  error  performance  were  given. 


Three  enhanced  DRDF  structures  were  presented.  Each  of 
these  new  approaches  used  parallel  filter  sections.  The  parallel 


-40- 


filters  could  share  much  of  the  basic  DRDF  hardware.  The  optimum 
approach  was  to  use  integer  value  taps  constructed  from  three 
parallel  sections.  The  need  for  hardware  multipliers  was  avoided 
through  the  use  of  shifting  and  adding. 

Expressions  for  the  RMS  error  of  the  integer  tap  method 
were  derived.  It  was  shown  that  the  RMS  error  is  inversely  pro- 
portional to  the  maximum  integer  used.  A low  pass  four  pole 
Chcbychev  filter  was  used  as  an  example.  Similar  results  have 
been  obtained  for  other  low  pass  structures  and  the  results  given 
in  this  paper  are  typical  of  what  might  be  expected. 
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Estimated 


Actual 


2.8  x 10 


1.9  x 10 


1.4  x 10 


1.1  x 10 


7.2  x 10 


5.7  x 10 


3.6  x 10 


3.0  x 10 


1.8  x 10 


1.5  x 10 


0.9  x 10 


0.8  x 10 


Comparison  of  Estimated  and  Actual 
Time  Domain  RMS  Errors 
for  the  Integer  Tap  Filter 


TABLE  II 
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Im 

Estimated 

Actual 

2 

0.396 

0.240 

4 

0.198 

0.140 

8 

0.102 

0.065 

16 

0.051 

0.033 

32 

0.025 

0.012 

64 

0.012 

0.007 

Comparison  of  Estimated  and  Actual 
Frequency  Domain  RMS  Errors 
for  the  Integer  Tap  Filter 


TABLE  III 


FIG.  2-  TYPICAL  LOW  PASS  FINITE  IMPULSE  RESPONSE 


AMPLITUDE 


FIG.  6-  IDEAL  MAGNITUDE  RESPONSE  AND  MAGNITUDE  RESPONSE  ERROR 
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Appendix  B 

A DIGITAL  FILTER  DESIGN  PROGRAM 
UTILIZING  THE  BILINEAR  Z TRANSFORM 

by 

Harriette  Markos 
and 

Thomas  A.  Brubaker 

Department  of  Electrical  Engineering 
Colorado  State  University 
Fort  Collins,  CO  80521 


ABSTRACT 

This  report  describes  the  design  of  digital  filters 
obtained  by  the  Bi^inear-Z  transformation  of  analog  filters. 
The  ideas  formulated  in  this  report  have  been  coded  into 
a computer  program  which  calculates  the  digital  transfer 
function  coefficients.  The  magnitude  function  is  then 
calculated.  The  Bilinear-Z  transform  is  applied  with  the 
frequency  warping  carefully  accounted  for  to  obtain  a 
realizable,  stable  digital  filter.  The  dependence  on  the 
sample  time  is  shown  in  the  comparison  of  the  digital  and 
analog  magnitude  functions.  The  phase  functions  are  not 
considered. 


This  work  was  supported  by  Contract  <fF33615-75-C-1138,  Air 
Force  Avionics  Laboratory,  Wright  Patterson  Air  Force  Base, 
Ohio  45433 
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I.  INTRODUCTION: 

In  network  analysis  the  transfer  function  is  a fundamental  character- 
istic of  a system.  The  transfer  function  H(s)  is  defined  by  McGillem  and 
Cooper  [1J  to  be  the  ratio  of  the  Laplace  Transforms  of  the  response  or 
output  signals  to  the  excitation  or  input  signals  when  the  initial  condi- 
tions are  zero.  If  the  excitation  voltage  is  represented  as  v^(t)  and 
the  response  voltage  is  represented  by  vQ(t),  then  the  transfer  function 


is  given  by 


H(s)  = 


L[vo(t)]  Vo(s) 
Ljvj(t) 1 = VjTsT 


(1) 


This  H(s)  transfer  function  is  of  interest  because  it  describes  the 
network  behavior.  For  simplicity,  it  will  be  referred  to  as  the  transfer 
function . 

This  report  and  the  computer  program  it  describes  are  concerned  with 
second  order  transfer  functions 


H(s) 


C s + C, s + C~ 
o 1 2 

2 

D s + D.s  + D,, 
o 12 


(2) 


These  second  order  structures  are  fundamental  building  blocks  for  lowpass 
filters,  highpass  filters,  bandpass  filters,  and  bandstop  filters.  For 
more  information  on  these  structures  see  Budak  [2]. 

Design  of  these  filters  is  often  done  by  specifying  the  magnitude 
function  H(w).  In  a filter  design  procedure  the  specifications  on  the 
magnitude  function  are  usually  concerned  with  critical  frequencies  such 
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as  j>dc.  the  frequency  at  DC;  and  the  frequency  where  the  response 

has  decreased  by  3db  from  the  DC  value.  For  high  Q filters,  the  fre- 
quency where  the  magnitude  function  peaks  is  typically  specified.  This 
frequency  is  referred  to  as  the  center  frequency,  wc.  These  frequencies 
are  the  critical  frequencies  under  consideration  in  this  report.  This 
report  is  written  to  outline  the  procedure  for  use  of  the  Bilinear-Z 
transform  as  applied  to  second  order  filters.  The  BLZ  program  is  a 
computer  program  to  perform  the  procedure  outlined  in  the  report. 

The  program  accepts  a second  order  function  in  s.  Before  apply- 
ing the  extended  Bilinear-Z  transform,  the  critical  frequencies  must 
be  prewarped  as  explained  in  Rabiner  and  Gold  [3].  Then,  the  extended 
Bilinear-Z  transform  described  as 

_ 2 Z-L  (3) 

S T Z+l 


is  implemented  and  the  equivalent  discrete  time  transfer  function  in  Z 
is  obtained: 


H (Z) 


+ AjZ  + A2 
+ BjZ  + B2 


(4) 


The  magnitude-squared  function  for  (4)  is  obtained  by  setting  Z » 
and  taking  the  sum  of  the  squares  of  the  real  and  imaginary  parts.  In 
the  design  program,  plots  of  both  the  analog  and  digital  filter  mag- 
nitude functions  are  available  for  presentation  on  a display  screen. 

The  program  recognizes  real  or  complex  poles  as  well  as  real  or 
complex  zeros.  When  equation  (2)  contains  both  zeros  and  poles,  the 
program  will  consider  two  cascaded  sections,  one  of  zeros  and  one  of 


poles 
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r 


H(s) 


D_  C s2  + C.s  + C„  C. 

2 . _o 1 2 . 2 

r D 2 

Z 2 D s + D.s  + D_ 

o 12 


where  the  cascaded  sections  may  be  considered  as 


hzro<s)  ■ 


C s + C. s + C- 
o 12 


and 


“P0L<S>  ' 


D s + D.s  + D„ 
o 12 


K HZRO^S^  ' ^POL 
(5) 

(5a) 


(5b) 


The  program  is  developed  with  specific  consideration  for  second  order, 
lowpass  filters.  However,  any  second  order  realizable  structure  can  be 
processed.  The  program  will  proceed  as  described  with  consideration  of  the 
three  critical  frequencies:  DC  frequency  ~3db  frequency  and  center 

frequency  u>c . Proper  operation  of  the  program  requires  the  filters  to 

have  an  or  an  ui  in  order  to  calculate  the  discrete  time  transfer 

3 c 

function  (4),  with  the  magnitudes  at  the  critical  frequencies  preserved. 

The  magnitudes  at  any  other  frequency  may  not  correspond  due  to  the  frequency 
warping. 

II.  DEVELOPMENT : 

A.  BILINEAR-Z  TRANSFORM 

The  extended  Bilinear-Z  transformation  is  described  in  Rabiner  and 
Gold  [3]  and  McGillem  and  Cooper  [4]  as  a mapping  from  the  s-plane  to  the 
Z-plane  with  the  following  properties:  the  s«jw  axis  is  mapped  onto  the 
unit  circle  of  the  Z-plane;  the  left  half  of  the  s-plane  (s<o)  is  mapped 
to  the  Interior  ol  the  unit  circle  In  the  Z-plane;  and  the  right  half  of 
the  s-plane  (s>o)  is  mapped  to  the  exterior  of  the  unit  circle  in  the 
Z-plane.  Thus,  the  analog  transfer  function  H(s)  given  by  (2)  is  trans- 
formed by  (3)  to  the  corresponding  digital  transfer  function  H(Z)  given 
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by  (4).  The  digital  transfer  function  H(Z),  evaluated  at  Z-e^U^,  is 


periodic  in  w with  period  oj 


2u 


P T 

The  extended  Bilinear-Z  transform  is  given  as 


, as  explained  in  McGillem  and  Cooper  [4] 


= 1 Z-l 
S T Z+l 


(3) 


and  is  applied  by  direct  substitution  for  s.  When  the  jw  axis  is  mapped 

onto  the  unit  circle  in  the  Z-plane,  there  is  a nonlinear  relationship 

between  the  analog  frequency  w and  the  corresponding  digital  frequency  52. 

Rabiner  and  Gold  [3]  illustrate  this  nonlinearity  as  follows: 

1QT 

Using  equation  (3)  let  s=jw  and  Z=e  so 

• - 2 ejOT-l  , 

T jflT 

eJ  +1 


then  multiplying  numerator  and  denominator  of  the  right  side  by  e 


-jOT/2 


jw 


2 ejOT/2_e-jOT/2 
j ejOT/2+e-jOT/2 


Recalling  the  exponential  relations  from  Euler's  equation 


2 . TAN  52T 
T 3 2 


and  so 


“ = T 


TAN  OT 
2 


(6) 


Therefore,  the  analog  frequency  m is  mapped  to  the  digital  frequency  51 
by  the  relationship  given  in  equation  (6).  When  it  is  desirable  to  have 
the  discrete  transfer  function  possess  the  same  magnitude  as  the  analog 
transfer  function  at  a particular  frequency  uj  D*  thC  analog  frequency 
imis t he  prewurped  or  altered  ho  that  when  the  Billnear-Z  transformation 
Is  perlormed  this  prewarped  frequency  will  map  into  the  desired 
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digital  frequency  The  is  determined  by 

2 TAN  ^D1  . (7) 

T 2 

The  relationship  of  this  to  the  coefficients  of  (4)  is  determined 
and  the  prewarped  coefficients  for  (4)  are  evaluated  using  the  prewarped 
value  for  ui^.  The  process  is  described  in  more  detail  in  the  following 

sections. 

B.  COMPLEX  CONJUGATE  ROOTS 

A second  order  structure  with  complex  conjugate  roots  is  considered 

as 


Ms  + M,s  + M « (s-r  )(s-r0) 
o 1 2 12 


(8) 


where  the  leading  coefficient  of  s is  normalized  to  1 by  dividing  through 
by  Mq.  The  roots  are  complex  conjugates 


r1  = a + jB 
r2  = a - j8  . 


(9a) 

(9b) 


Equation  (8)  may  also  be  considered  as 


, A ,2  ^ 2 2 o . 2 

(s+a)  + 8 = s + q—  s + a)o 


(10) 


where  iu  - 


•V? 


and 


rt 

o . U) 

Q'^  ° 


In  equation  (10),  represents  the  distance  from  the  origin  to  the  location 
of  the  roots,  and  Q indicates  the  slopes  of  the  radial  lines  that  connect 
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the  root  locations  to  the  origin.  For  a<<8,  Q>>1. 


£ 


Figure  1 


a,(5,u)o,  Q are  related  as  follows: 


a = 


2Q 


f i - — 2 
6 = a >o\/  4Q^ 


J7-  - 


+ B 


sfJ-  - 


+ 6 


2a 


A plot  of  the  magnitude  function  H(w)  for  complex  conjugate  poles 
reveals  a peaking  effect  with  the  maximum  value  at  the  center  frequency, 


<*>c.  This  is  covered  extensively  by  Budak  [5].  The  amount  of  peaking  is 


indicated  by  Q,  where  Q may  be  considered  as  a measure  of  the  ratio 


DC^magnitude^6  * when  Q is  greater  than  five  the  peak  will  become  prom- 


inent so  that  the  center  frequency  wc  may  become  more  important  than 


the  -3db  frequency  These  high  Q lowpass  filters  actually  have  a 

bandpass  type  structure.  Therefore,  it  is  Important  to  preserve  the 
appropriate  critical  frequency.  The  program  utilizes  a user  input  limit- 
ing value,  Q-limit.  If  Q is  less  than  or  equal  to  Q-limit,  the  -3db 
frequency  is  prewarped.  But  if  Q is  greater  than  Q-limit,  the  center 
frequency  u>j  is  prewarped. 


1 
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Expanding  and  simplifying: 


To  solve  for  w , apply  the  quadratic  formula  employing  the  positive  root 
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At  this  time  the  digital  frequency  can  be  prewarped  by  equation  (6) 
to  obtain  the  prewarped  analog  -3db  frequency  10^  where 

„ _ 2 TAN  ^ . 

0)3  - x 2 


Now  utilizing  equation  (11)  and  w ',  the  prewarped  coefficients  for 
H(s)  can  be  evaluated.  Recognizing  first  that 


D.  = w 
2 o 


n ^5 

1 " Q " Q 


then  using  (12) 


2(u)')2 

°2  (2  - i )2  + n/(2  - -,)2  +4 

Q Q 


D'  = 1 . 
o 

C2  C2 

To  evaluate  C~,  note  the  DC  gain  at  H(s“0)  is  DC  * — - “ —j 

z o 


Since  it  is  critical  to  retain  the  DC  gain,  C'  is  evaluated  to  be 


dc  • d; 


The  prewarped  coefficients  determined  by  (15),  (16),  (17),  and  (19) 
now  form  a new  transfer  function  H'(s)  which  when  transformed  by  the 
Bilinear-Z  transformation  will  exactly  match  in  magnitude  the  original 
transfer  function  H(s)  at  the  critical  frequencies  ci)^,  anc*  . 
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Instead  of  prewarping  co^,  the  center  frequency  wc  could  be  prewarped. 
Again,  it  is  first  necessary  to  locate  a>c . Recalling  that  ioc  is  the 
center  frequency  of  the  peaking  effect,  differential  calculus  is  used 
to  find  the  location  of  this  extrema.  Taking  the  derivative  with  respect 
to  to  of  the  magnitude-squared  function  for  (11) 

H(ui=toc) 


d 

dto[  | 


I (<o  2-<o  2)2  + (“o“c) 

L ° c ~ 
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Now,  employing  equation  (6)  the  prewarped  analog  center  frequency 

uj"  is  obtained: 

c 


U) 

c 


2 TAN 
T 2 


The  prewarped  coefficients  are  found  from  (13),  (14),  (17),  (19), 
and  (20)  to  be: 


(21) 


(16) 


D'  = 1 
o 


(17) 


C'  = DC  • 


°2 


(19) 


As  in  the  case  of  prewarping  w y these  prewarped  coefficients 
determine  a new  transfer  function  H'(s)  which  when  transformed  by  the 
Bilinear-Z  transformation  will  exactly  match  in  magnitude  the  original 
transfer  function  H(s)  at  the  critical  frequencies  and  u> c . 

Example  1 : 

Consider  H(s)  * — r , where  Q « 10. 

s + 0.1s  + 1 

Using  a sample  time  of  T * 1.0, 

u>  = 0.997497  and  = 1.551026. 
c 3 

The  digital  transfer  function  coefficients  can  be  obtained  on  print  out 
during  execution  of  the  BLZ  program. 
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Figures  2 and  3 are  plots  of  the  analog  and  digital  magnitude 
functions  with  prewarped  and  u>c  prewarped,  respectively.  In  both 
cases,  a sample  time  of  1.0  second  was  used.  The  sample  time  affects 
the  H(Z)  magnitude  function:  For  small  T values  (less  than  0.1 
second)  the  two  magnitude  functions  are  nearly  the  same  for  all  u> 
values.  As  the  sample  time  increases,  the  two  magnitude  functions 
begin  to  differ  from  each  other  for  all  w values  except  the  critical 
frequencies.  The  sample  time  of  1.0  second  shows  this  variation 
clearly.  Values  of  T greater  than  1.0  second  start  to  show  larger 
differences  to  the  point  that  distortion  causes  the  two  magnitude 
functions  to  differ  in  an  unacceptable  manner. 

2.  Complex  Conjugate  Zeros 

Recalling  equations  (2),  (5a),  and  (10)  the  transfer  function 
with  complex  conjugate  zeros  can  be  considered  as 


1 

rs2 

+ clS  + C2' 

1 

r 2 
s 

ID  0 

+ s + ID  3 

Q o 

°2 

1 J 

' °2 

_ 

1 J 

where  the  s coefficient  is  . irmalized  to  1 by  division  throughout  by 

C . 
o 

Considerations  for  the  -3db  frequency  and  the  center  frequency 
o)c  follow  immediately  from  the  development  for  complex  conjugate  poles 
and  are  not  reiterated  here.  The  effect  of  Q in  this  case  is  to  cause 
a dip  in  the  magnitude  function  and  the  u>c  is  found  by  locating  a min- 
imum rather  than  a maximum. 

The  prewarped  coefficients  for  H(s)  are  found  as  follows: 

I.  Prewarped  to 

"3  ’Mr 

2 TAN  w3T 

“ t 2 


(2  - 1 ) + 1(2  - l )2  + A 

Q V Q 


(22) 
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2^) 

(23) 

(2  - ^2>  + ^<2  " ^2>2  + * 

[cT 

\ 2 

(24) 

VQ 

« 1 

(25) 

c: 

2 

(26) 

II.  Prewarped  wc: 


(27) 

(28) 

(24) 

(25) 

(26) 


In  both  cases  the  prewarped  coefficients  form  a new  transfer  function 

H'(s)  which  when  transformed  by  the  Bilinear-Z  transformation  will  exactly 

match  the  original  transfer  function  H(s)  at  the  critical  frequencies 

io > and  w-  in  case  I;  or  w and  <*>  in  case  II. 

DC*  j DC*  C 

C.  REAL  ROOTS 

A second  order  structure  with  real  roots  is  considered  as 

Ms2+M,8+M  « (s-r.Ms-r,) 

o 1 2 1 i 


(27) 
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where  the  roots  are  evaluated  by  application  of  the  Quadratic  formula 


r.  ■ -M  +7m2  - AM  M. 
1 1 v 1 o 2 

2M 


(28) 


-M, 


^7 


4M 

o 2 


(29) 


2M 


When  complex  conjugate  roots  are  under  consideration  it  is  first 
necessary  to  locate  the  critical  or  wc  frequencies.  In  the  real 
root  case  this  is  not  necessary.  If  each  root  is  prewarped  and  the 
prewarped  H(s)  coefficients  determined  from  the  prewarped  roots,  the 
new  transfer  function  H'(s)  will  yield  an  H(z)  after  the  Bilinear-Z 
transform  which  matches  in  magnitude  the  original  H(s)  at  the  DC 
frequency  to^  and  at  one  other  frequency  dependent  on  the  sample  time 
T used.  The  technique  is  to  factor  the  filter  coefficients  to  locate 
the  roots.  Then  each  root  is  prewarped  and  the  sections  multiplied 
together  before  applying  the  Bilinear-Z  transform. 

The  roots  are  found  by  equations  (28)  and  (29) . Then  using  equation 
(6),  the  prewarped  roots  are 


r T 

2 TAN  1 


r T 

2 TAN  _2_ 
r2  = T 2 


(30) 

(31) 


Then  the  prewarped  coefficients  are  found  by  (27),  (28),  and  (29) 


to  be 


M'  = 1 
o 


(32) 


M'  - -(^  + r2) 


(33) 
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M2  ’ rir2  * 

(34) 

If  M1 
0 

0 there  is  only  one  real  root 

-Mo 

2 

r’^ 

2 TAN  rT 

(35) 

Then 

r = T 2 

and  so 

M"  = 0 

0 

(36) 

M£  = 1 

(37) 

h:  = -r ' • 

(38) 

A special  case  arises  when  M2  = 0 and  M^  “ 0,  or  when  * 0 and 
M2  = 0,  causing  two  roots  at  zero,  or  one  root  at  zero,  respectively. 
In  both  cases  TAN  (0)  = 0 indicating  a linear  relationship,  hence,  no 
prewarping  is  required. 

1.  Real  Poles 

Referring  to  equations  (2)  and  (5)  the  transfer  functions  with 
real  poles  can  be  considered  as 


H(s) 


1 

+ D^s  + D2 


(39) 


where 


D •*-*+! 
o o 


1 

[2 


in  the  previous  development. 
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The  roots  (poles)  and  prewarped  coefficients  can  be  evaluated  as 
described  before  with  C'  being  determined  from 


hence 


C'  = DC  • D'  . 


2.  Real  Zeros 

Referring  to  equations  (2)  and  (5)  the  transfer  function  with  real 
zeros  can  be  considered  as: 


H(s>  -i  jV+ClS  + C2] 


where 


C +-**! 
o o 


Cl~*l 


C2-M2 


in  the  previous  development. 


The  roots  (zeros)  and  prewarped  coefficients  can  be  evaluated  as 

C2 

described  with  D'  defined  by  (40)  to  be  ^ . (41) 


Example  2 : 


Consider 


s + 15s  + 50 


with  poles  at  5.0  and  10.0. 


Figure  4 is  a plot  of  the  analog  and  digital  magnitude  functions  for 
this  lowpass,  real  pole  filter  using  two  sample  times.  For  the  sample 
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time  of  0.3  sec.,  the  analog  and  digital  magnitude  functions  agree 
exactly  at  uip^.,  and  at  an  intermediate  frequency  However,  for 

the  sample  time  of  0.2  sec.,  the  analog  and  digital  magnitude  functions 
agree  at  w , and  at  a frequency  ui„  different  from  u>  . When  evaluating 
real  root  structures,  the  sample  time  has  a direct  effect  on  the 
critical  frequency  value  because  the  roots  are  prewarped,  not  the 
critical  frequency.  Note  that  for  first  order  sections  prewarping 
the  root  is  equivalent  to  prewarping  the  -3db  frequency.  For  the 
first  order  sections  in  cascade  the  net  -3db  frequency  will  not  match 
because  it  was  not  prewarped. 

D.  H(Z)  TRANSFER  FUNCTION 

After  the  prewarping  is  accomplished,  the  digital  transfer  function 
is  evaluated.  The  transformation  is  made  by  direct  substitution  of  (3) 
into  the  transfer  function  of  prewarped  coefficients  giving 


Multiplying  numerator  and  denominator  by  [T(Z+1)1 

C ' [ 2 (Z-l) ] 2 + c:[2(Z-l)][T(Z+l)]  + C'[T(z+l)]2 
__  v o 1 o 

H(Z)  = F'TT(^T)1TT  d;[2(z-i)][t(z+i)]  + d'[t(z+i)]2 

o 1 o 


Expanding,  and  then  collecting  terms 

H(Z>  ■ |4Co  * 2TCl'  * t2c2'z2  * I2t2c2'  - 8cq1Z  * - 2TC{  + t2c2>  (41) 


[ 4D ' + 2TD,'  + T2D']Z2  + [2T2d: 

OIL  L 


8D']Z  + [4D'  - 2TD.'  + t2d;] 
o o 1 2 


Recalling  equation  (4) 


AQZ  + A^Z  + A2 


Z + + b2 


H (Z) 


(4) 
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where 


and 


DEN  ■=  4D'  + 2TD.'  + T2d: 

o 1 2 


A = [4C'  + 2TC:  + t2c:]/den 

0 o 1 l 

A,  = [2T2C'  - 8C']/DEN 

1 l O 

A.  = [4C'  - 2TC;  + t2c;]/den 

l O 1 2. 

B,  - [2T2d;  - 8D']/DEN 

1 L o 

B„  = [4D'  - 2TD;  + t2d:]/den 

z o 1 l 


e.  calculation  of  magnitude  functions 

1.  Analog  Transfer  Function  H(s) 

The  magnitude  function  for  the  analog  transfer  function  is  obtained 
by  setting  s=jto  in  (2)  and  finding  the  square  root  of  the  sum  of  the 
squares  of  the  imaginary  and  real  parts: 


Hs(w) 


| H (s-jw) 


Cq  (jto)  + Cx  (j<o)  + C2 


Dq  (Jo>)  + Dx  (jto)  + D2 


and  since  j * \j-l 

(C2-CoU)2)  + j(C1u>) 


H (to) 
s 


(D2-Dq(0  ) + j (Djto) 


Then  the  magnitude  of  a quotient  is  the  quotient  of  the  magnitudes 


or  equivalently 


. / (C.-C  (o2)2  + (C.io) 

«»  - Vi-  2 °22  -1- 

W(VV>  + 


(D^to) 


2 


(C  -C  <o2)2  + (C.  to) 
H >>  -J  2 ° 1 

(U2"Do'1’  } + (0i,j) 


(42) 


2.  Digital  Transfer  Function  H(Z) 

The  magnitude  function  for  the  digital  transfer  function  requires 
special  consideration  of  the  original  analog  transfer  function.  If  H(s) 


4 


is  composed  of  one  section  of  either  poles  or  zeros,  the  digital  magnitude 
function  is  found  using  the  prewarped  H'(s)  transfer  function.  But,  if 
H(s)  is  composed  of  two  sections  in  cascade  as  in  equation  (5),  there  are 
two  transfer  functions  to  be  considered:  the  transfer  function  of  pre- 
warped coefficients  for  the  pole  section,  H^Q^(s);  and  the  transfer 
function  of  prewarped  coefficients  for  the  zero  section,  The 

resultant  total  transfer  function  is  the  product  H^^Cs)  * Hzro^’ 

The  technique  used  here  is  to  transform  each  section  to  obtain  the  digital 
transfer  function  in  cascades  H'  (Z)  and  H'  (Z) , then  evaluate  the 

lUL  L KU 

digital  magnitude  functions  for  each  section,  and  the  resultant  total 
digital  magnitude  function  is  the  product  Rp0L(w)  • HZR0^*  In  both 
cases  the  general  technique  to  obtain  the  magnitude  function  H(u>)  from 
a digital  transfer  function  H(Z)  is  to  let  Z * e^uT,  apply  Euler's 
formula,  and  find  the  square  root  of  the  sum  of  the  squares  of  the 


imaginary  and  real  parts 


Hz(u)  = |H(Z  = ejuT) 


Recalling  Euler's  formula 


A ej2wT  + A eju)T  + A, 

;t2  ' 


eJ  = cos  wT  + j sin  wT 

and  combining  real  and  imaginary  parts  yields 

[A  cos  2u)T  + A, cos  u)T  + A.  ] + 1 [A  sin  2u>T  + A,  sin  wT ] 
Z^'  ~ [ cos  2(jjT  + B1cos  coT  + B2]  + j[  sin  2u>T  + B sin  uT] 


and  since  the  quotient  of  the  magnitudes  is  the  magnitude  of  the  quotients 


V“> 


[A  cos  2uiT  + A,  cos  u>T  + A_]2  + [A  sin  2u)T  + A,  sin  u>T]2 
o 1 / o 1 

2 2 
[ cos  2wT  + B^cos  idT  + B21  + [ sin  2uT  + B^sin  u>T] 


In  the  case  of  the  analog  transfer  function  being  composed  of  two 


sections  in  cascade,  the  resultant  total  digital  magnitude  function  is 
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found  as  Che  product  of  the  magnitude  functions  of  each  section.  This 
multiplication  effectively  destroys  the  correspondence  of  the  critical 
3db  or  center  frequencies  for  each  section,  resulting  instead  in  a new 
critical  frequency.  Only  if  the  critical  frequency  of  prewarping  is 
identical  in  both  sections  will  the  resultant  critical  frequency  be  the 
same.  However,  if  critical  frequencies  for  the  total  transfer  function 
are  prewarped,  such  as  in  Butterworth  filter  design,  then  the  individual 
sections  are  not  treated  separately  in  terms  of  prewarping.  This  allows 
the  critical  frequencies  to  be  handled  directly  from  input  to  output. 
This  form  of  design  has  been  treated  in  two  separate  reports  now  being 
revised.  These  are  "A  Fortran  IV  Design  Program  for  Butterworth  and 
Chebychev  Band-Pass  and  Band-Stop  Digital  Filters",  and  "A  Fortran  IV 
Design  Program  for  Low  Pass  Butterworth  and  Chebychev  Digital  Filters". 
The  revisions  will  be  completed  in  early  August,  1976. 

III.  OPERATION  OF  THE  PROGRAM 

The  principles  developed  so  far  are  utilized  in  the  BLZ  program. 
Written  in  DEC  Fortran  IV,  this  program  accepts  a second  order  analog 
transfer  function,  performs  the  necessary  prewarping  to  match  critical 
frequencies,  implements  a Bilinear  Z transform,  and  determines  the 
digital  transfer  function  by  the  methods  explained  in  section  II. 

A.  INITIALIZATION 

The  operation  of  the  BLZ  program  assumes  a high  degree  of  user 
interaction  via  an  external  teletype  or  keyboard  device.  The  user  is 
asked  to  supply,  through  input,  the  following  initialization  factors: 

1.  The  program  will  request  input  specifying  real  poles  or 
> complex  poles,  real  zeros  or  complex  zeros.  Input  is  in 

the  form  of  1 for  yes  and  0 for  no.  As  the  program 
receives  each  response  a flag  is  set  to  each  affirmative 
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reply.  If  flags  are  set  for  both  zero  and  pole  sections 
a special  flag  is  set  to  indicate  two  sections  in  cascade. 

2.  The  program  will  request  the  range  of  w values  the  mag- 
nitude functions  are  to  be  plotted  over.  A minimum  m 
and  a maximum  u>  are  input  in  213  format.  This  option 
facilitates  the  opportunity  of  investigating  any  specific 
region  of  the  magnitude  function. 

3.  The  program  will  then  request  the  total  number  of  points 
to  be  plotted.  Input  is  in  13  format,  however,  only 
values  up  to  500  are  allowed  under  the  current  dimension 
allotments.  If  more  than  500  points  are  desired,  the 
dimension  statement  for  the  arrays  used  must  be  altered 
to  a value  greater  than  or  equal  to  the  total  number  of 
points  plus  two. 

4.  The  program  will  request  a sample  time  T in  seconds  and 
the  coefficients  of  the  analog  transfer  function  H(s). 

As  stated  before: 


H(s) 


C s2  + C.s  + C. 
o I 2 

D s2  + D.s  + D„ 
o 1 L 


(2) 


The  format  statement  here  allows  values  up  to  F9.5.  The 
leading  coefficients  need  not  be  normalized  to  1. 

NOTE : All  input  values  must  be  presented  in  the  format  specified  or 
with  decimal  points  and  commas  included  in  all  relevant  positions. 


B.  EXECUTION 

After  the  initialization  factors  are  received,  the  program  proceeds 


to  calculate  the  analog  magnitude  function  defined  by  equation  (42). 
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The  user  may  have  the  computer  output  the  DC  value  before  responding  to 
a request  to  display  a plot  of  Hg(<D)  versus  u)  to  be  generated  on  a 
display  screen.  The  IDIOT  subroutine  is  called  for  this  purpose,  and 
is  included  in  the  program  package  as  explained  under  the  section  on 
software. 

The  program  control  now  advances  into  the  transformational  section 
of  its  operation.  If  the  special  flag  to  indicate  two  sections  in 
cascade  is  set,  the  program  will  evaluate  a magnitude  function  Hr0L(<d) 
for  the  pole  section,  then  a magnitude  function  H^rqOd)  for  the  zero 
section.  The  resultant  total  digital  magnitude  function  H (id)  is  the 

Lt 

product  of  and  HZrq(id).  If  the  special  flag  is  not  set,  the 

program  checks  first  for  the  pole  flags  and  then  for  zero  flags.  In 
all  cases,  when  a set  flag  is  encountered,  the  control  branches  to  the 
appropriate  subroutine  for  the  set  flag.  The  subroutines  serve  to 
evaluate  the  prewarped  analog  transfer  function  coefficients  as  described 
in  section  II.  In  the  case  of  complex  conjugate  poles  or  complex  con- 
jugate zeros,  the  user  is  requested  to  input  a limit  value  for  Q to 
determine  prewarping  of  tD^  or  idc-  In  all  cases,  the  user  will  have  the 
option  to  print  out  such  information  as  the  -3db  frequency,  the  center 
frequency,  or  the  prewarped  roots. 

As  each  prewarping  subroutine  is  completed,  control  returns  to  the 
main  program.  The  Bilinear-Z  transform  is  executed  as  described  and  the 
digital  transfer  function  is  obtained.  The  magnitude  function  defined 
in  (43)  is  evaluated  for  the  <d  range  indicated  by  the  initialization. 

A plot  of  this  magnitude  function  is  available,  and  if  the  user  so 
indicates,  the  digital  magnitude  function  will  then  be  superimposed  on 
the  plot  of  the  analog  transfer  function.  If  the  user  so  desires  the 
program  will  now  output  the  minimum  and  maximum  values  of  the  magnitude 
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f unction. 

At  this  time  the  user  is  afforded  the  option  to  alter  the  sample 
time  T.  If  he  responds  affirmatively,  the  program  requests  a new  T 
value  to  be  input.  The  control  then  branches  to  Inquire  if  the  user 
desires  to  replot  the  analog  magnitude  function.  If  the  user  responds 
negatively,  the  present  display  screen  contents  are  retained  and  the 
next  plot  will  be  superimposed  with  these  same  plots.  If  the  user 
responds  affirmatively,  the  screen  is  cleared  and  the  analog  magnitude 
function  is  displayed  and  the  program  proceeds  as  before. 

In  the  case  of  two  sections  in  cascade  the  program  evaluates  for 
each  section,  then  offers  a plot  of  the  resultant  total  digital  magnitude 
function.  The  options  for  the  minimum  and  maximum  magnitude  values  and 
a new  sample  time  T are  offered  as  described  before. 

When  the  user  is  satisfied  with  the  sample  time  used  and  rejects  the 
option  to  alter  T,  the  program  continues  with  an  offer  to  printout  the 
table  of  values  of  the  magnitude  functions  over  the  range  of  u values. 

At  this  point,  in  the  case  of  two  sections  in  cascade,  the  user  is 
offered  the  opportunity  to  display  plots  of  the  magnitudes  of  each  section 
individually. 

Two  more  printout  options  are  available  at  this  time:  the  prewarped 
normalized  analog  transfer  function  coefficients;  and  the  digital  transfer 
function  coefficients.  The  final  option  is  to  run  the  program  over  again. 
An  affirmative  response  directs  control  to  the  Initialization  portion  of 
the  program,  while  a negative  response  terminates  the  program. 

C.  HARDWARE 

The  BLZ  program  is  written  in  DOS  Fortran  version  9.02.  It  was 
developed  on  a PDP-11/20  with  a DOS/BATCH  operating  system.  The  function 


plots  were  obtained  using  a GT40  graphics  display  terminal,  with  hardcopy 
plots  available  on  a Houston  Instruments  Plotter  Interfaced  to  the  GT40. 
All  printed  results  were  available  using  a Centronics  line  printer.  The 
program  is  written  to  be  easily  modified  for  use  with  systems  of  similar 
conf iguration. 

D.  SOFTWARE 

Operation  of  the  BLZ  program  requires  considerable  user  interaction 

for  data  input.  The  program  writes  instructions  and  questions  to  unit  6. 

The  user  responds  with  the  appropriate  data  which  is  read  from  unit  4. 

These  two  units  must  be  assigned  to  appropriate  output  and  input  devices 

at  run  time.  Additional  output  is  available  to  unit  5,  which  must  also 

be  assigned  to  a device  at  run  time.  The  program  was  written  and  tested 

with  units  4 and  6 assigned  to  a teletype  keyboard  and  unit  5 assigned 

to  a line  printer.  A sample  run  using  these  devices  is  provided  in  figure 

The  BLZ  program  is  composed  of  a main  block  with  four  subroutines  to 

perform  the  prewarping: 

RZRO 

RPOL 

CZRO 

CPOL 

and  a subfunction  TAN  to  calculate  the  trigonometric  function  tangent. 

The  program  also  calls  the  following  routines  from  library  files: 

FLOAT  IDIOT 

SQRT  RANGE 

ABS 

COS 

SIN 

The  routines  on  the  left  are  called  from  the  DOS/BATCH  FTNLIB.  The 
routines  on  the  right  are  required  for  plotting.  They  are  called  from  the 
PLTLIB  file  and  are  included  in  the  program  package. 
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The  IDIOT  subroutine  was  written  for  plotting  applications  on  the  GT40 
graphics  display  terminal.  Any  subroutines  required  by  IDIOT  are  included 
in  PLTLIB,  along  with  the  GT40  plotting  routine  (PDP-11  assembly  language 
MACRO).  Explanation  and  documentation  of  the  use  of  IDIOT  can  be  obtained 
by  listing  the  program.  All  of  the  required  plotting  software  for  the  GT40 
has  been  included  with  the  program  package,  or  the  user  may  incorporate  his 
own  plotting  routines. 

Recommended  set  up  procedures  for  Use  of  the  BLZ  program  are  as  follows: 

1.  Compile  BLZ  program.  This  assures  system  compatibility. 

2.  Compile  PLTLIB. SRC. 

3.  Assemble  SENDGT.  This  is  the  plotting  routine  to  run  the 

GT40. 

4.  Build  the  subroutine  library  PLTLIB. OBJ  including 

1.  PLTLIB. OBJ 

2.  SENDGT. OBJ 


in  that  order. 

5.  Assemble  and  link  PLOTGT.MAC  for  latter  use  in  GT40  plotting. 

6.  Link  together  BLZ,  PLTLIB,  FTNLIB. 

7.  Load  PLOTGT.LDA  into  GT40  and  start  it  running. 

8.  Run  BLZ  (any  changes  on  I/O  devices  mav  be  made  now  by  assign- 
ment statements  for  the  devices  4,  5,  and  6 as  explained 
previously. ) 


Figure 
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INTRODUCTION 

This  report  describes  how  to  use  two  programs  for  the  weighted 
least  squares  design  of  nonrecurslve  and  recursive  digital  filters. 
First  the  theoretical  aspects  are  considered  and  the  design  equations 
are  developed.  The  signal  model  In  this  work  la  assumed  to  he  a 
polynomial  because  the  state  model  Is  simple.  However,  the  theory  is 
easily  extended  to  Include  any  signal  model  represented  by  a linear 
differential  equation. 

Then  the  operation  of  the  two  programs  is  described  along  with 
examples  to  illustrate  their  operation. 
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LEAST  SQUARES  DESIGN  OF  NONRECURSIVE  DIGITAL  FILTERS 

The  design  of  nonrecursive  and  recursive  digital  filters  using 
weighted  least  squares  is  based  on  a model  for  the  input  signal. 

The  most  common  models  that  are  currently  in  use  are  d l f t erent I a I 
equations.  The  subsequent  representation  of  the  differential  equa- 
tion by  a fir9t-order  vector  differential  equation  leads  to  the 
concept  of  a state  variable  and  the  state  space  representation  for 
a system.  By  use  of  the  proper  formulation,  a continuous  signal 
represented  by  a differential  equation  can  be  described  by  a first- 
order  vector  difference  equation  or  a discrete  time  state  space 
model.  References  that  describe  the  essential  aspects  of  state 
variables  are  by  De  Russo,  Roy  and  Close  [1]  and  Chen  ( 2 J . 

Since  signals  from  laboratory  Instruments  are  not  usually  des- 
cribed by  a differential  equation,  approximations  often  utilize  a 
polynomial.  For  this  reason,  the  vector  form  of  a polynomial  approx- 
imation will  be  used  in  this  report.  The  reader  should  be  aware 
that  this  can  be  generalized  to  include  any  signal  model  that  can  be 
represented  as  a linear  time-varying  differential  equation.  Later 
in  the  paper  a scalar  model  representing  a Gaussian  time  signal 
will  also  be  utilized  to  develop  a time-varying  filter  that  can  be 
used  for  reducing  the  base-line  error  and  for  the  initial  separat ion 
of  signal  components. 

To  develop  the  polynomial  model  let  the  signal  z(t)  be  repre- 
sented by  a polynomial  of  order  m.  At  the  time  t =*  nT  the  state 
vector  for  the  signal  is  given  by 


the  use  of  a Taylor  series  representation  for  each  element  of  x(nT) 
now  permits  the  state  of  system  at  t*(n+h)T  to  be  described  in  terms 
of  the  state  at  t=nT  by  the  relationship 

x[ (n  f h)T]  • ♦Ih |x[nTl  (1) 

where  <t>[h]  is  the  mxm  state  transition  matrix  with  elements 

*ijlhl  " i ! ( fey  hJ  l (4) 

i<J<m 

- 0 KJ 

The  state  transition  matrix  41 1 h ] satisfies  all  of  the  relationships 
for  general  state  transition  matrices  with  the  Important  two  being  for 
this  work 

-1 


♦ l— h ] * t ♦ (h)  ] 


(5) 


<t> [ m ] <t>[p]  * $lm  + pl  (b) 

At  this  point,  these  models  can  be  utilized  in  the  design  process. 

De sign  o f Nonrecurslve  Fllte rs 

Let  the  input  signal  start  at  time  t=*0  and  assume  that  the  signal 
over  a finite  data  window  is  approximated  by  a polynomial  z(t). 

Defining  the  state  of  the  signal  by  (2),  then  the  state  of  the  signal 
at  time  t=(n+h)T  is  given  in  terms  of  the  signal  at  time  t=nT  by  (3). 
Given  that  the  signal  starts  at  t=0,  the  first  S observations 
are  each  defined  as 

m(jT]  = Mx  [ j T 1 + y ( l T J j-0,1 ,2, . . . t.-l  (7) 

where  M is  a row  matrix  that  relates  the  measurable  state  variables 
to  the  actual  measurements.  The  elements  of  each  noise  vector  v[jT| 
are  the  measurement  noises.  In  general  these  are  taken  as  random 
variables  with  zero  mean  and  the  time  dependent  autocovariance  matrix 

R [ j T ] = E(v(jT)yt(jT) ] (8) 

However,  in  most  laboratory  systems  only  the  data  is  available  and 
the  derivatives  are  not  measurable.  Furthermore,  the  noise  covariance 
matrix  is  usually  not  known  and  for  scalar  measurements  the  noise 
variance  is  assumed  constant  and  the  noise  samples  uncorrelated.  For 
the  remainder  of  this  paper,  only  scalar  measurements  and  uncorrelated 

measurement  noise  with  time -invariant  statistics  will  be  assumed.  The 

2 

mean  value  of  the  noise  will  be  taken  as  zero  and  the  variance  as  o . 
The  results  are  easily  extended  to  vector  measurements  and  to  measure- 


ment  noise  with  time  varying  statistics. 
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For  i.  observations,  the  total  observation  vector  at  t»nT  is 


defined  as 


n>t  [ nT  1 


m[nT| 
ml (n- l )T ) 


M[ (n-f+l)T ] 

J 

This  vector  now  forms  a data  window  of  !.  data  points.  For  an  estimate 
of  the  data  at  t=nT  the  use  of  the  expression 


x[(n-j)T|  = $(-j)  x[nT] 


is  combined  with  (9)  to  yield 


M$(-l) 


mt[nT]  = 


x[nT  ) + v^ 1 nT] 


M4>[-y.-K  ] 


The  matrix  of  constants  H[nT]  is  now  defined  as 


H[nT]  = 


M4>  (-1 ) 


M4>[-)t+l] 


so  that  (11)  can  be  written  as 


m.  [nT]  = H[nT]  x[nT]  + v [nT] 


The  elements  in  H[nT]  are  constants  given  by 


(-i)j  0<i<l 


(14) 
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where  l is  the  size  of  the  data  window  and  q is  the  order  of  the 
polynomial  plus  one.  When  i=j«0  the  value  of  (14)  is  one.  Note  that 
since  (12)  is  a matrix  of  constants  there  is  no  need  to  make  the 
matrix  a function  of  time.  However,  in  the  derivation  for  the  recur- 
sive filters  this  provides  a method  for  separating  different  H matrices. 

The  optimal  estimate  of  the  data  at  t=nT  is  now  given  in  terms 
of  weighted  least  squares  or  minimum  variance  because  this  form  is 
utilized  in  the  derivation  of  the  recursive  filters.  If  the  covari- 
ance matrix  for  the  total  observation  vector  is 

Rt(nT)  = E[vt  (nT)  v^(nT)]  (15) 

the  optimal  minimum  variance  estimate  is 


x[nT]  = W[nT ] mt[nT] 

where  W[nT]  is  a series  of  constant  weights  given  by 


W[nT] 


Ht(nT)[Rt(nT)]~1H(nT) 


"1Ht(nT) [ R t (nT) )_1 


(16) 


(17) 


For  uncorrelated  noise  with  constant  variance  a (15)  is  a diagonal 
matrix  with  elements  a*  and  (17)  reduces  to 


W[nT ] = [Ht(nv)  H(nt)]-1  Ht(nT) 


(18) 


This  is  the  same  result  obtained  using  conventional  least  squares  when 
the  noise  covariance  is  a diagonal  matrix  of  equal  constants.  The 
reader  should  be  aware  that  the  estimate  vector  Is  .in  estimate  of  the 
data  and  all  of  the  derivatives  in  the  model.  However,  all  of  the 
weights  in  the  derivative  terms  must  be  scaled  by  the  scale  factors  in 
the  state  vector  defined  by  (2). 


1 


I 


-101- 


The  covariance  matrix  for  the  estimate  given  by  (16)  is  given  by 

S(nT)  = W(nT)  R^  (nT)  Wt(nT)  (19) 

Substituting  (17)  into  (19)  now  yields 

S (nT)  = [ H*1  (nT)  [R^nT)]"1  H(nT)]_1  (20) 

which  simplifies  further  to 

S(nT)  = [Hc (nT)  H(nT)]-1  a2  (21) 

for  the  uncorrelated  noise. 

By  defining  a delay  or  prediction  factor  a,  estimates  of  the  data 
otT  units  behind  or  ahead  of  the  point  t=nT  can  be  done.  To  do  this 
the  total  observation  vector  is  written  as 


m (nT)  = H (nT)  x[(n-a)T]  + v (nT) 

— t Ot  “ T 


(22) 


where  H (nT)  now  takes  the  form 
a 


H (nT)  = 

a 


M$(-n) 

M$(-a-l) 


1 


M<t>(-H.-a+l) 


(23) 


The  individual  elements  of  the  matrix  now  become 


Ha(nT)  i j = (a-i)j  0<i<*. 

0<j<q 


(24) 


The  form  for  the  optimal  estimate  now  utilizes  H^(nT)  in  (18).  The 

covariance  of  the  estimate  uses  H (nT)  in  (21). 

a 
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Example 

For  a five-point  window  with  a=0,  the  optimal  weight  matrix  and 
the  covariance  matrix  for  the  estimate  are  shown  in  Table  1 for  a 
third-order  \ jlynomial  fit.  For  a=2,  the  estimate  of  the  data  in 
the  middle  of  the  window  is  obtained  and  weight  matrix  and  covariance 
matrices  are  shown  in  Table  2.  This  case  corresponds  to  weights 
given  in  reference  [3]. 

In  practice,  the  design  of  nonrecursive  filters  is  initiated  by 
specifying  the  size  of  the  window  which  is  the  number  of  rows  in  the 
H matrix,  the  order  of  the  polynomial  approximation  which  is  one  less 
than  the  number  of  columns  in  the  H matrix  and  a which  determines 
the  coefficient  values  in  the  H matrix.  This  makes  the  design  suit- 
able for  use  with  interactive  graphics  since  only  three  parameters 
need  be  specified  to  generate  the  weight  and  covariance  matrices. 

If  the  model  of  the  signal  process  is  modified  by  additive 
uncorrelated  driving  noise,  the  variance  terms  of  the  driving  noise 
fade  the  memory  so  that  past  data  has  less  effect  on  the  estimate. 

This  can  be  thought  of  as  uncertainty  in  the  signal  model.  The  con- 
cept is  particularly  important  in  recursive  filter  design.  For  non- 
recursive filters,  it  can  also  be  utilized  and  can  be  useful  when 
using  a non-recursive  filter  to  initialize  a recursive  filter. 

If  the  model  includes  driving  noise,  the  state  at  time  t=nT  is 
given  by 

x[nT ] = M - 1 1 xl  (n-l)T  1 + w[(n-l)T]  (?■>) 

where  w[(n-l)T]  is  a sample  from  a noise  process  with  moan  zero  and 
covariance  matrix  Q.  This  noise  process  is  assumed  to  be  white.  Ln 


J 
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terms  of  the  total  measurement  vector,  m^tnT]  now  becomes 

m [nT]  = H[nT]  x[nT]  + n(nT)  + v [nT] 

— t - *■%  — t 

where  £t[nT]  Is  the  total  noise  vector  due  to  the  driving  noise, 
scalar  measurements  this  is  given  by 

r° 

| -M<J> [ -1  ] w[  (n-l)T ] 

j^tnT]  = ; -M4>  f-2  ] wf  (n-1  )T  ] - M4>(-1]  w[(n-2)Tj 

I 

I • 

The  total  noise  vector  that  corrupts  the  total  measurement  vector  is 
now  defined  as 

r t [ nT  ] = |>t  [nT  ] + yt[nT]  (28) 

For  scalar  measurements  the  covariance  matrix  for  rt[nT]  has  diagonal 
elements  whose  value  increases  down  the  diagonal.  Since  the  diagonal 
elements  are  not  the  same,  the  minimum  variance  expressions  of  (16)  and 
(17)  must  be  employed  to  find  optimal  linear  estimates. 

Example 

To  illustrate  how  the  driving  noise  affects  the  filter  weights, 
consider  a zero  order  process  which  is  equivalent  to  estimating  a signal 
of  constant  value.  For  a three-point  filter  with  scalar  measurements, 
the  total  noise  vector  is 

r 

j v[nTj 

v[(n-])T|  - w[ (n-1 )T J 
v[(n-2)T]  - w[ (n-2 )T ] - v[(n-l)T] 


(26) 

For 


(27) 


r [nT] 


(29) 
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If  the  variance  of  the  measurement  noise  is  and  of  the  driving 

a 

ise  the  covariance  matrix  of  rt(nT]  is 


no 


Rt[nT] 


o 0 0 

0 a2  + a,2  0 
0 0 


2 2 
o + 2o^ 


The  resulting  weight  matrix  obtained  using  minimum  variance  is 


where 


W[nT]  * [aQ  a^  a2l 


, 2 2.  2 2. 

(o  + a ) (2o,  + o ) 

; v — ' 

I i 


*1 


and 


where  D is  given  by 


D 


2 2 2 
0(0^  + o ) 


D - (2o^2  + a2)  (0l2  + o2)  + a2  (2 a2  + o2) 


+ a2(a^2  + o2) 


(30) 


(31) 

(32) 
(3!) 

(34) 

(35) 


The  terms  in  the  weight  matrix  satisfy  the  following  inequality 

a0  i al  a2 

For  o^2  “ 0,  the  equalities  hold  and  ag  “ ai  * a2  = which  are 

2 

the  well  known  weights  to  estimate  a mean.  If  o^  is  not  zero,  the 

2 2 
inequalities  hold  and  as  becomes  larg?>  with  respect  to  o , a^ 

approaches  one  and  a^  and  a2  become  smaller  so  that  fading  is 


(36) 
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introduced.  The  covariance  of  the  estimate  Is  a scalar  given  by 

. (a.2  + (2a.  ^ + a2)  o2 

S(nT)  « ^ - Co2  (37) 

where  C is  a dimensionless  constant  whose  value  approaches  one  as  o^2 

2 2 2 

becomes  larger  than  o . Thus  for  >>  o only  the  current  observa- 

tion Is  effectively  used  and  the  variance  of  the  estimate  Is  approximately 
2 

o . 

The  fading  obviously  reduces  the  signal-to-noise  enhancement  of  a 
nonrecursive  digital  filter.  On  the  other  hand,  fading  can  be  used  to 
reduce  the  deterministic  error  due  to  a nonexact  model.  in  practice, 
fading  in  nonrecursive  filters  is  not  often  utilized.  Instead,  the 
window  length  is  more  typically  used  as  a design  parameter.  In  future 
work,  however,  it  may  be  desirable  to  further  explore  the  relationship 
between  fading  and  the  size  of  the  data  window  to  achieve  improved  designs 
when  nonexact  signal  models  are  employed. 

LEAST  SQUARES  DESIGN  OF  RECURSIVE  DIGITAL  FILTERS 
The  fixed  memory  or  nonrecursive  filter  design  is  now  extended  to 
recursive  digital  filters  that  utilize  all  of  the  data.  The  result  is 
a recursive  form  that  is  usually  called  the  Kalman  filter.  While  there 
are  a variety  of  derivations  for  the  Kalman  filter,  starting  with  a fixed 
memory  filter  using  a polynomial  model  with  driving  noise  gives  the  result 
in  such  a way  to  give  the  reader  greater  intuition  about  how  the  filter 
works. 

The  derivation  of  the  recursive  filter  is  started  by  using  a signal 

model 

x[nT  ] - 4>  [ 1 ] x[  (n-l)T  ] +w[(n-l)T) 


(38) 
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and  the  scalar  observation  or  measurement  model 

y [nT]  = Mx[nTl  + v[nTj. 


( 19) 


In  (38)  and  (19)  the  terms  v|(n-l)T|  .ind  v | n I | /ire  I lie  » I r I v 1 1 1 k imWc 
and  the  measurement  noise.  These  noise  terms  have  zero  mean  and  ire 
uncorrelated  with  themselves  and  each  other.  Given  the  n+1  measurements 
starting  at  t=0,  the  total  observation  vector  at  time  t=nT  is  given  in 
terms  of  x(nT]  as 

y[nT]  = H[nT]  + £t[nT]  + VjJnT]  (40) 

In  (40)  the  matrix  H[nT]  is 


I n'l'  I 


| - 1 | 


M<J>(-n] 


the  vector  n [nT]  is 


*1 


— M<t>  [ -1  ] w[  (n-1  )T  ] 


-M.^  4>[-(3-j)  ] w [ (n-j)T] 


-Mj^  <t' [ - (n+l-j  ) ] w[(n-j)T] 


and  v^JnT]  is  the  total  measurement  noise  vector.  Defining  the  sum 


(41) 


(42) 


r t [ r.T  ] = pt[nT]  + yt[nT] 


(41) 
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as  the  total  noise  vector  wilh  autocovar lance 


R Am]  = E{r  InT]  r InTl) 
t “ t t 

the  optimal  estimate  using  linear  minimum  variance  is  given  by 


x[nT]  = W[nT]  ^[nT] 


where  W[nT]  is  the  weight  matrix 


W[nT] 


H*"  [nT^R^  [nT]  _1  H[nT]  ,_1  H* [nT]  [r^ [nT ] _l 

L J - L 


The  covariance  of  the  estimate  is 


S[nT)  >■  W[nT]  RtlnT]  W [nT] 
Substitution  of  (46)  and  (47)  gives  an  alternate  form 


S [nT]  - HL[nT]  Rfc[nT]  H[nT] 


i 

H[nT]j 


which  allows  the  alternate  form  for  (46)  to  be  written  as 


W[nT]  = S [nT]  Ht[nT]  Rt  [nT]  1 


The  forms  given  by  (46),  (47),  (48)  and  (49)  apply  to  the  remaining 
estimates  used  in  the  derivation  and  the  reader  should  remember  them. 

Next  the  prediction  or  forecast  of  the  signal  state  X[(n+1)T]  is 
found  by  first  writing  the  total  observation  vector  ^[nT]  as 


Yt[nT]  = H1[nT]  x[(n+l)T]  + pu[nT]  + v^nT]  (50) 


where  H^[nT]  is  given  as 


H^nT]  - H[nT]  4>l -1 J 


and  £^[nT]  is 
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-M<t(-1]  w[nT  1 

2 

-ME  f[-(3-j)]  w|(n-l-))T) 


%tlnl]  * 


n+1 

-M  E 4>[-(n+l-j  ) 1 w[  (n-l-j  )T 

j=l 


(52) 


The  term  y^JnT]  is  the  same  as  in  (40)  since  no  more  observations  have 
been  taken. 

The  estimate  of  x[(n+l)T]  is  now  given  by 


XjKn+lJTl-  WjfnTj  £t[nT] 


(53) 


where 


W1fnT]  = 


H^nT] 


RuInT) 


-1 


H^nT] 


r n _i 

H tnT]|R , JnT]! 

L J 

(54) 


In  (54)  is  the  covariance  matrix  of  the  nclse  vector 


rlt(nTl  = Elt  fnT J + v^InT) 


(55) 


The  corresponding  covariance  matrix  for  the  estimate  x^((n+l)T]  is 


S1[(n+1)T]  = W1[nT]  R^lnT]  wJ(nT] 


H^lnTl 


Rlt(nT] 


-1  H-i 

H.[nT] 

1 J 


(56) 


The  development  of  the  relat ionship  between  R [nTJ  and  Rt(nT]  is 
the  next  issue. 

To  obtain  this  recognize  that  since  the  random  variables  w|nT) 

and  v[nT]  are  uncorrelated,  R^fnTJ  is  the  sum  of  the  covariance 

matrices  for  p^fnTj  and  vfc[nT].  For  the  latter,  if  the  scalar  noise 

2 

v(nT]  has  variance  o , the  covariance  matrix  of  yt[nT]  is  a diagonal 
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2 

n+l  x n+1  matrix  with  elements  o . Taking  the  covariance  of  £^t(nT] 
and  performing  some  algebraic  manipulation  gives 

E{£u(nT]  fnT] } = Ef^lnT]  £[[nT]}  + H^nT]  QH^  [nT ] (57) 

where  Q is  the  diagonal  covariance  matrix  of  the  driving  noise  vector 
w[nT],  This  matrix  is  usually  assumed  to  be  time  invariant.  Thus,  from 
(57) 

Rlt[nT]  = Rt[nT]  + H^nT]  QH^  [nT  ] (58) 

Substituting  (58)  into  (56)  now  gives 

S1((n+1)TJ  = W^nTj  R^nTl  wJ[nT]  (59) 

+ W1[nTj  H1[nT]  QH* [nT ] wJ[nT] 

If  X^[(n+1)T]  is  an  unbiased  estimate  then  the  constraint  relationship 

W^nT]  H^nT]  = I (60) 

must  be  satisfied  [4]  so  that  (59)  can  be  simplified  to 

SjKn+DT]  - W [nT]  Rt[nT]  wJ[nT]  + Q (61) 

Also  recognizing  that  x^[(n+l)T]  can  be  written  as 

XjKn+DT]  = *[1]  x[nT ] (62) 

the  weight  matrix  W^[nT]  is 

W^nT]  = 4>[1 ) W(nT) . (63) 

Substitution  at  (58)  into  (61)  and  applying  (47)  now  gives  the  final 
form  for  S^[(n+1)T]  as 

S1f(n+l)Tl  =•  | 1 J S[nT]  <tt[nT]  + Q (64) 

The  recursion  is  now  formulated.  When  the  observation  at  t * (n+l)T 
arrives  the  new  total  observation  vector  in  terms  of  the  signal  state 


Mi 
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x[(n+l)T]  is 

y 1 1 ( n-*- 1 ) T ] - H|(n+1)T|  + pt|(n+l)T|  + vtl(n+l)T|.  (6r>) 
The  estimate  la  given  by 


x[(n+l)T]  - S[(n+1)T]  [ (n+l)T ] Rjfn+DTl'1 


with  the  covariance 

S[(n+1)T] 


Ht[(n+1)T]  Jr^ [ (n+l)T]]  _1  H[(n+1)T] 


(66) 


(67) 


where  Rt[(n+1)T]  is  the  covariance  matrix  of  the  total  measurement  noise 


vector 


rt[(n+l)T]  = pt[(n+l)T]  + y^Un+DT] 


(68) 


First  the  recursion  for  the  covariance  matrix  is  found.  Given  that 


-M4>[ -1 1 w[nT] 


£tt(n+l)T] 


I -M  !:  [ — ( 3—  i ) 1 w[  (n+l-j  )T  ] 


n+1 

-M  l 4> [ — (n+L— j ) ] w[  (n+l-j  )T ] 
j-1 


(69) 


substitution  of  (52)  into  (69)  gives 

0 
I 

£ [(n+l)T] 


2ltlnT] 


(70) 


Next  H[(n+1)T]  is  given  as 


H[(n+1)T] 


M 

M<b  [ -1  ] 

^ M<t>[-(n+l)  | 


(71) 
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whlch  can  be  rewritten  with  the  aid  of  (41)  as 


H | nT  | - <t>l-l 

H[nT| 


The  covariance  matrix  for  v^[(n+l)T)  is  now  an  n+2  x n+2  diagonal 

2 

matrix  with  elements  a so  that  Rt[(n+1)T]  is  seen  to  be 


o2|  0 


Rt[(n+1)T]  = -- 


0 | Ru[nT] 


Since  this  is  a diagonal  matrix  its  inverse  is 


2 i 0 


Rt[(n+1)T] 


i ^ r i -i  ! 

0 i jvnTj  i 


Substitution  of  (74)  and  (72)  into  (67)  and  performing  the  matrix 

multiplication  now  yields 

— ~] 

S{(n+1)T]  + <t>C [ -1 1 HC[nT]  R1(.[nT]  _1  H[nT|  4>( -1  ] _1  (75) 

a L 

Substitution  of  (51)  into  (75)  and  the  use  of  the  inverse  of  (56)  gives 
~ M^M  —1  ” ^ 

S[(n+1)T]  - ~ + Sl(nT]  (76) 

o I — — I 


Equation  (76)  along  with  (64)  now  forms  a recursion  for  the  covariance 
of  the  estimate. 
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Kor  the  quantities  In  (76)  the  application  of  the  matrix  Inversion 
lemma  [5],  gives 


M^M 

2 + 

Sl(nTl 

-1 

0 

-- 

_ 

-l 


SjlnTl 


(77) 


- ;;  [nT|  MC[M  S j | nT  | Ml  + n‘|'  M S^nTl 


This  form  will  be  used  to  generate  an  expression  for  the  Kalman  gain.  Post 
multiplying  both  sides  of  (75)  by  MC/a2  and  using  (75)  the  Kalman  gain  is 
defined  as 


K[(n+1)T]  - S[(n+1)T]  Mt/o2 


(78) 


S.[(n+1)T]  M^o2  + M S^tn+DT]  MC  ]_1 


Thus  the  covariance  matrix  given  by  (76)  becomes 


S[(n+1)T]  *=  [I  - K[  (n+1  )T  ] M|  S1l(n+l)Tl 


(79) 


The  recursion  for  the  optima)  estimate  Is  now  formed  that  uses  the  Kalman 
gain  given  by  (78).  Using  the  form  for  the  optimal  estimate  given  by 
(47)  the  estimate  x[(n+l)T]  is 


-1 


Y.t  l (n+l)'r ] 


x[(n+l)T]  = S [(rt-l)T ] Ht[(n+l)T]:RIt[(n+l)Tl 

I 

where  ^t[(n+l)T]  is  the  new  total  observation  vector  given  by 


(80) 


£t((n+l)T]  = 


y[(n+l)T? 


(81) 


Zt  lnTl 

Substitution  of  (74),  (72)  and  (81)  into  (80)  and  carrying  out  the  matrix 
multiplication  yields 


x[(n+l)T]  - S[ (n+l)T ] 


y[ (n+l)Tl 

2 

o 


(82) 


+ <t>  ( -1  J H [ nT  | 


V"T| 


■l 


ytlnT) 
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Using  the  estimate  of  the  forecast  x^((n+l)T)  (82)  is  simplified  to 


x[(n+l)T]  - S[(n+1)T] 


~ y[(n+l)T!  + I Sx[ (n+l)T) 


-1 


XjJ  (n+1 )T ] 


Adding  and  subtracting  (MtM/o'^)  x^|(n+l)T],  performing  some  algebraic 
manipulation,  and  applying  (76)  now  yields 


x[(n+l)T)  = xx[ (n+l)T]  + K[(n+1)T] 


(83) 


y[(n+l)T]  - M x^  [ (n+l)T  ]j  (84) 

where  K[(n+1)T]  is  given  by  (78).  Equations  (84),  (79),  (78),  (64)  and 
(62)  now  form  the  recursive  formulation  called  the  Kalman  filter.  These 
equations  are  now  summarized  as 


x^nT)  = 4>(1)  x[  (n-l)T ] , (85) 

S^nT)  = $(1)  S [ (n-l)T]  ^(l)  + Q , (86) 

K(nT)  * S1(nT)  M^o2  + M S^nT)  M^-1  , (87) 

S(nT)  - [f  - K(nT)  Mj  S^nT)  (88) 

and 

x(nT)  = x^(nT)  + K(nT)  [y(nT)  - M x^nT)]  (89) 


- [I  - K(nT)M]  4>(1)  x[  (n-l)T ] + K(nT)  y(nT) 

In  this  set  of  equations,  x^(nT)  is  the  forecast  or  prediction  of  the 

estimate  at  t«nT  using  the  previously  generated  estimate  at  t«|(n-l)T| 

The  covariance  of  the  forecast  is  S^(nT).  The  term  K(nT)  Is  the  t Imo 

varying  Kalman  gain  matrix  and  y(nT)  is  the  observation  at  t«nT.  The 

terms  x(nT)  and  S(nT)  are  the  estimate  at  t»nT  and  its  covariance. 

2 

The  term  a is  the  variance  of  the  measurement  noise  and  the  term  Q Is 
the  covariance  of  the  driving  noise.  It  is  the  term  Q that  serves  as 


a key  design  parameter. 
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If  there  is  no  driving  noise  so  that  the  diagonal  elements  ol  1} 
are  all  zero,  the  filter  Is  simply  an  expanding  memory  filter.  This  means 
that  if  the  filter  is  initialized  properly,  the  estimates  will  correspond 
to  those  obtained  by  designing  nonrecursive  filters  where  the  data  window 
starts  at  zero  and  a new  weight  matrix  W(nT)  is  computed  as  each  new 
measurement  is  made.  Obviously  as  the  window  expands,  the  variance  of 
the  estimate  will  decrease;  however,  if  the  model  is  not  exact  determinis- 
tic errors  will  begin  to  Increase. 

In  using  the  equations,  they  must  he  initialized  properly  if  a trulv 
unbiased  estimate  is  to  be  formed.  In  practice  this  is  usually  done  using 
a nonrecursive  filter.  For  exact  initialization,  driving  noise  must  be 
included  in  the  computation  of  the  nonrecursive  filter  weights.  To  mini- 
mize the  computation,  the  minimum  H matrix  should  be  used  which  means  the 
number  of  rows  should  equal  the  number  of  columns.  By  using  this  minimum 
H matrix  the  filter  weights  are  computed  and  the  initial  optimal  estimate 
is  formed  from  the  actual  data. 

USING  THE  PROGRAMS 

HARDWARE 

The  programs  are  written  in  DOS  FORTRAN.  They  were  developed  on  a 
PDP  11/20  with  a DOS/BATCH  operating  system.  Printed  results  are  written 
to  logical  unit  5,  which  can  be  assigned  at  run  time  to  a line  printer, 

CRT  terminal,  disk  file,  or  other  suitable  output  device.  Data  is  entered 
from  units  6 and  3,  which  can  be  assigned  to  a card  reader,  disk  data 
file,  TTY  keyboard,  or  any  other  suitable  input  device.  Plots  can  he 
obtained  with  a GT40  graphics  display  terminal  and  the  plotting  subroutines 
provided.  The  programs  can  be  easily  modified  to  use  other  plotting  routines. 
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NONRECURSIVE  FILTER  PROGRAM 


The  program  for  generating  the  coefficients  for  non-recursive  filters 
is  called  PROGRAM  WINDOW.  WTNDOW  is  written  in  DEC  FORTRAN,  but  may  be 
run  on  other  versions  of  FORTRAN  IV  with  minor  modifications.  The  program 
can  call  plotting  packages  to  produce  CRT  or  hardcopy  plots  of  the  filter 
response  to  various  inputs. 

Program  WINDOW  reads  the  filter  parameters  SIGMA,  N,  M,  IA,  IPLOT 
from  logical  unit  6,  where: 

SIGMA  determines  the  input  covariance  matrix  R. 

If  SIGMA  > 0.0,  R becomes  o2*I,  where  o2 
= SIGMA  and  I is  the  identity  matrix. 

If  SIGMA  < 0.0,  the  matrix  R is  read  from 
unit  6 in  10F8.2  FORMAT. 

N is  the  number  of  points  in  the  window 

(1  < N < 20). 

M is  the  order  of  the  polynomial  fit  (1  < M < 9) . 

IA  is  the  offset  a from  the  first  point  in  the 

window.  If  IA  > 0,  the  filter  predicts  IA 
sample  times  ahead  of  the  most  recent  sample. 

If  IA  < 0,  the  filter  smooths  IA  sample  times 
behind  the  most  recent  sample. 

IPLOT  is  the  plotting  control  variable.  Tf  IPLOT 

* 0,  the  program  finishes  after  the  coeffi- 
cient matrices  are  computed  and  printed 
out.  If  IPLOT  0,  an  input  signal  is  read, 
and  the  input  points  are  filtered  using  the 
coefficients  computed. 

The  parameters  are  read  in  F10.4,  413  FORMAT. 

Once  the  filter  parameters  have  been  read,  WTNDOW  generates  the  re- 
quired S (covariance)  and  T matrices  , and  computes  the  weight  matrix  W. 
All  three  matricies  are  then  written  to  Logical  unit  5.  If  IPLOT  » 0. 
the  program  terminates  after  the  weight  matrix  W has  been  printed. 

If  IPLOT  j 0,  WINDOW  reads  101  sample  points  from  logical  unit  3 in 
G15.6  FORMAT.  These  points  are  provided  by  the  user,  and  are  used  by 
WINDOW  as  the  filter  input  signal.  The  program  filters  this  input  signal 
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using  the  weight  matrix  W,  and  tabulates  the  input  signal,  output  signal 
(filter  response),  and  error  signal  (Input  minus  output)  for  each  sample 
time.  The  tabulated  results  are  printed  on  unit  5.  The  sum  of  the  total 
absolute  error  is  also  computed  and  written  to  unit  5.  WINDOW  then  calls 
the  plotting  package  routines  (subroutine  IDIOT)  to  plot  the  input  and 
output  (estimate)  signals  vs.  time.  After  a PAUSE,  the  program  calls  the 
plotting  package  to  plot  the  error  signal  (input  minus  output)  vs.  time. 
The  plotting  packages  are  included  to  be  used  with  WINDOW,  or  WINDOW  mav 
be  changed  to  call  other  plotting  routines. 

Program  WINDOW  can  be  modified  to  write  the  derivative  estimates. 
Note  from  equation  2 that  the  m1"*1  derivative  term  is  multiplied  by  a 
Tm/m!  factor,  where  T is  the  sample  time.  Thus,  if  the  derivative 
estimates  are  written  out,  they  will  be  scaled  by  this  factor.  To 
illustrate  the  use  of  the  program,  WINDOW  was  run  with  the  parameters 
SIGMA  = 1.0,  N = 5,  M = 3,  IA  = 0,  and  IPLOT  = 0.  The  filter  weighting 
coefficients  are  listed  row  by  row,  and  are  shown  with  the  input  para- 
meters and  S and  T matrices  in  Fig.  1.  The  filter  obtained  using  the 
weight  matrix  shown  estimates  the  input  data  and  the  first  three  deri- 
vatives. In  equation  form,  these  estimates  are  given  by 


x[nT]  = 0.9857  y[nT]  + 0.05714  y[(n-l)Tj  - 0.0857  y!(n-2)T] 
+ 0.05714  y[(n-3)T]  - 0.01429  y[(n-4)T]  , 

x[nT  | - 1.488  y[nT]  - 1.619  ^LLn-iJTJ.  “ -r>/l4  yl(n-2_)T| 

r 

+ 1.048  y [ (n-  i)T J - 0.3452  y[(n-4)T] 

T ’ 

x[nT]  - 0.6429  y[nT]  - 1.071  y[(n-l)T]  - 0.1429  y[(n-2)T] 

t2/2 

+ 0.9286  y [ (n-3)T ] - 0.3571  y[(n-4)T) 


-117- 


x[nT]  =*  0.08333  y[nT]  - 0.1667  y[(n-l)T] 


+ 0.1667  y[(a-3)T]  - 0.08333  y[(n-4)T] 


T3/6 


RECURSIVE  FILTERS 

The  Kalman  filter  program  Is  called  ADAPT.  ADAPT  Is  written  In  DEC 
FORTRAN,  but  may  be  run  on  other  versions  of  FORTRAN  IV  with  minor  modi- 
fications. The  program  can  call  plotting  packages  to  produce  CRT  or 
hardcopy  plots  of  the  filter  response  to  various  inputs. 

Program  ADAPT  reads  the  filter  parameters  SIGMA,  M,  and  Q from 
logical  unit  6,  where 

SIGMA  determines  the  initial  covariance 

matrix  R.  If  SIGMA  > 0.0,  R becomes 
o^*l,  where  * SICMA  and  I is  the 
identity  matrix.  If  SIGMA  < 0.0, 
the  matrix  R is  read  from  unit  6 in 
10F8.2  FORMAT. 

M is  the  order  of  the  polynomial  fit 

(1  < M < 9). 

Q is  the  driving  noise  term. 

The  parameters  are  read  in  F10.4, I3,F10.4  FORMAT. 

Once  the  filter  parameters  have  been  read,  ADAPT  generates  an  initial 
S (covariance)  and  T matrix.  ADAPT  then  generates  an  initial  weight 
matrix  W,  which  is  used  to  initialize  the  x vector.  The  initial  covari- 
ance matrix  is  used  to  initialize  the  S matrix. 

Once  initialized,  ADAPT  reads  101  sample  points  from  logical  unit  3 
in  G15.6  FORMAT.  These  points  are  provided  by  the  user,  and  are  used  by 
ADAPT  as  the  Kalman  filter  input.  The  program  filters  the  input  using 
equations  85  through  89.  The  sample  number  filter  input,  filter  output, 
error  signal  (input  minus  output),  and  Kalman  gain  are  tabulated  and 
written  to  unit  5.  The  total  absolute  error  is  also  computed  and  written 


A 
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to  unit  5.  ADAPT  then  calls  the  plotting  package  routines  (subroutine 
IDIOT)  to  plot  the  Input  and  output  (estimate)  signals  vs.  time.  After 
a PAUSE,  the  program  calls  the  plotting  package  to  plot  the  error  sig- 
nal (input  minus  output)  vs.  time.  After  a second  PAUSE,  IDIOT  is 
called  to  plot  the  Kalman  gain  vs.  time. 

Program  ADAPT  can  be  modified  to  write  the  derivative  estimates. 
Note  from  equation  2 that  the  m^1  derivative  term  is  multiplied  by  a 
Tm/m!  factor,  where  T is  the  sample  time.  Thus,  if  the  derivative 
estimates  are  written  out,  they  will  be  scaled  by  this  factor. 

To  illustrate  the  use  of  the  program,  ADAPT  was  run  with  the 
parameters  SIGMA  = i.O,  M = 3,  Q = 0.0.  The  S (covariance)  and  T 
matrices  used  to  generate  the  initial  weight  matrix  (W)  are  shown  in 
Fig.  2.  The  S and  W matrices  are  used  to  initialize  the  Kalman  fil- 
ters' S and  x matrices.  The  Kalman  filter  was  then  used  to  filter 
an  ideal  sinusoidal  signal.  A portion  of  the  tabulated  results  is 
shown  in  Fig.  3. 

GETTING  ON  LINE 

In  order  to  run  program  WINDOW  and  ADAPT,  first  build  two  files 
named  WINDOW. FTN  and  ADAPT. FTN  from  the  sources  provided  (card  deck 
or  paper  tape).  Also  create  PLT.FTN  and  SENDGT.MAC  from  the  sources. 
Compile  programs  WINDOW  and  ADAPT  with  the  FORTRAN  compiler  to  create 
WINDOW. OBJ  and  ADAPT. OBJ.  The  one  word  integer  option  should  be 
selected  for  all  compilations.  Also,  compile  PLT.FTN  to  create 
PLT.OBJ.  Assemble  SENDCT.MAG  under  the  MACRO  assembler  to  create 
SENDGT.OBJ.  Create  a subroutine  library  called  Pl.TLIB.OBJ  from 


PLT.OBJ  and  SENDGT.OBJ  (in  that  order).  Next,  LINK  WINDOW. OBJ, 
PLTLIB.OBJ  and  the  FTN  library  to  create  a file  called  WINDOW. LDA. 


Also,  create  ADAPT. LDA  by  LINKing  ADAPT. OBJ,  PLTLIB.OBJ  and  the  FTN 
library.  The  files  WINDOW. OBJ,  ADAPT. OBJ,  PLT.OBJ  and  SENDGT .OBJ  mav 
now  be  deleted. 

Before  running  WINDOW  or  ADAPT,  build  the  file  PLOTCT. MAC  from 
the  source.  PLOTCT  is  the  plotting  routine  that  is  loaded  into  the 
CT40.  Assemble  and  LINK  PLOTCT  so  that  PLOTCT. LDA  may  be  loaded  into 
the  GT40.  After  PLOTGT.LDA  is  running  in  the  GT40,  ADAPT. LDA  or 
WINDOW. LDA  may  be  executed  using  a RUN  command.  The  source  listings 
contain  additional  documentation  on  these  programs. 
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***FINITE  MEMORY  DIGITAL  FILTER  PACKAGE ♦ ♦ ♦ 


VARIANCE  OF  ERROR'  1 0000 
SIZF  OF  WINDOW-  3 
ORDER  OF  FIT-  2 

TO  PREDICT  0 UNITS  AWAY  FROM  THE  FRONT  OF  THE  WINDOW 


MAGIC  T MATRIX 


1.  0000 

0 0000 

0.  0000 

0 0000 

1 0000 

-1.  0000 

1 0000 

-1  0000 

1.  0000 

-2.  000 

4.  000 

-8  000 

1 0090 

-2  000 

9 000 

-27.  00 

1 000O 

-4  000 

16.  00 

-64  00 

VARIANCE  MATRIX,  S 


0.  9357 

1.  488 

0.  6429 

0 

3222E-01 

1 488 

6.  279 

2.  869 

0 

5972 

0 6429 

2 869 

2.  571 

0 

4167 

0 3222E-01 

0 3972 

0 4167 

0 

6944E-01 

WINDOW  WEIGHTS, 

FROM  FRONT  TO  BACK 

0 9857 

ROW  1 OF 
0.  3714E-01 

W,  CORRESPONDING  TO  THE  INPUT  SIGNAL 
-0  S371E-01  0 5714E-01  -0  1429E-01 

1 488 

ROW  2 OF 
-1.  619 

W,  CORRESPONDING  TO  DERIVATIVE  NUMBER 
-0  5714  1 043  -0  2452 

0 6429 

ROW  2 OF 
-1.  071 

W,  CORRESPONDING  TO  DERIVATIVE  NUMBER 
-0  1429  0 9236  -0.  2371 

ROW  4 OF  W,  CORRESPONDING  TO  DERIVATIVE  NUMBER  2 
0 3222E-01  -0  ±€€7  -0.  2576E-06  0 1667  -0  S222E-01 


Fig.  1 - Program  WINDOW  Sample  Output 
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***KALMAN  STRUCTURE  DIGITAL  FILTER  PACKAGE*** 

DRIVING  NOISE3  0 0000 
ORDER  OF  FIT-  2 

INITIAL  ERROR  VARIANCE  1 0000 


1 0000 
1 0000 
1 0OO0 
1 0000 


MAGIC  T MATRIX 
0.  0000 
-1  0000 
-2.  000 
-2  000 


0 0000 
i 0000 
4.  000 
9.  000 


0 0000 
-1  0000 
-8  000 
-27.  00 


VAR I ANCE  MA  TR I X,  S 


1 0000 

1 822 

1 0000 

0 1667 

1 822 

14  72 

12.  50 

2.  611 

i 0000 

12  50 

11  50 

2.  500 

0 1667 

2 611 

2.  500 

0.  5555 

WINDOW  WEIGHTS, 


FROM  FRONT  TO  EACH 


1 0000 


ROW  1 OF  W,  CORRESPONDING  TO  THE  INPUT  SIGNAL 
0 4128E-05  -0  2502E-05  0 9527E-06 


1 


ROW  2 OF  W, 
-2.  000 


CORRESPONDING  TO  DERIVATIVE  NUMBER 
1 500  -0  2222 


1 


1 0OO0 


0 1687 


ROW  2 OF  W, 
—2.  500 


CORRESPONDING  TO  DERIVATIVE 
2 000  -0  5000 


ROW  4 OF  W.  CORRESPONDING  TO  DERIVATIVE 
— 0 5000  0 5000  —0.  1667 


NUMBER 


NUMBER 


Fig.  2 - Program  ADAPT  Sample  Output 


TINE 

OUTPUT 

ESTIMATE 

ERROR- 

iRL 

GH ! N 

4 

000 

0 

2SS4 

0 2394 

0 

26S2E-06 

0 

9857 

5. 

000 

0 

4794 

0 4794 

0 

1627E - 04 

0 

9695 

»£. 

000 

0 

5646 

0 5646 

0 

5460E-04 

0 

^ ;2;  O 

~7 

000 

0 

6442 

0 6441 

0 

1017E-02 

0 

2162 

O 

000 

0 

7174 

0 7172 

0 

1563E-07 

0 

3724 

9 

000 

0 

7822 

0 7821 

0 

2229E-02 

0 

3419 

10 

000 

0 

3415 

0 3412 

0 

2122E-07 

0 

3057 

11 

00 

0 

9912 

0.  3908 

0 

4220E-02 

0 

7707 

12 

00 

0 

9220 

0 9214 

0 

6002E-02 

0 

7282 

1.2 

00 

0 

9626 

0 9627 

6i 

S295E-02 

0 

7079 

14 

00 

0 

9855 

0 9842 

0 

1140E-02 

0 

6797 

15 

00 

0 

9975 

0 9959 

0 

1551E-02 

0 

6524 

16 

00 

0 

9996 

0 9975 

0 

2086E-Ci2 

0 

6239 

1? 

00 

0. 

9917 

0 9389 

0 

276SE-02 

111 

6061 

IS 

00 

0 

9728 

0 9702 

0. 

2620E-02 

0 

5843 

19 

00 

0. 

9462 

0.  9416 

0 

4669E-02 

0 

5649 

20 

001 

0 

9092 

0.  9024 

0 

5923E-02 

0 

5462 

21 

00 

0 

3622 

0 8556: 

0 

7 450E— 02 

0 

5287 

V 

00 

0 

8085 

0.  7992 

0 

9229E-02 

0 

5122 

22 

00 

0 

7457 

0 7244 

0 

1120E-01 

0 

4 9 6 8 

24 

00 

0. 

6755 

0 6613 

0 

1267E-01 

0 

4322 

25 

00 

0 

5935 

0 5821 

0 

1625E-01 

0 

4634 

26 

00 

0 

5155 

0 4961 

0 

1927E-01 

0 

4552 

9"? 

00 

0. 

4274 

0.  4047 

0 

2272E-01 

0 

4420 

00 

0 

2250 

0 20:36 

0 

2640E-01 

0 

4212 

y q 

00 

0 

2292 

0 2088 

0 

2040E-01 

0 

4201 

"0 

00 

0 

1411 

0 1064 

0 

2472E-01 

0 

4095 

21 

00 

0 

4158E-01 

0 2265E-02 

0 

2922E-01 

0 

2925 

22 

00 

-0 

5S27E-01 

-0  1025 

0 

4417E-01 

0 

TOOCl 

-•  -■ 

00 

-0 

1577 

-0  2070 

0 

4924E-01 

6i 

2807 

24 

00 

-0 

2555 

-0  2100 

0 

5447E-01 

61 

2720 

25 

00 

-0 

2503 

-0  4106 

0 

5981E-01 

0 

2627 

26 

00 

-0 

4425 

-0  5077 

0 

6519E-01 

0 

2557 

00 

-0. 

5293 

-0  6004 

0 

7054E-01 

0 

2480 

00 

-0. 

6119 

-0.  6876 

0 

7576E-01 

6i 

2407 

2’ 9 

00 

-0 

6373 

-0  7685 

0. 

3076E-01 

0 

T;  "7 

40 

00 

-0 

7563 

-0  3422 

0 

3542E-01 

0. 

2269 

41 

00 

-0. 

3132 

-0  9080 

0 

3968E-01 

0 

2204 

• • 

00 

-0 

3716 

-0  9649 

0 

9226E-01 

6i 

2142 

4" 

00 

-0 

9162 

-1  012 

0 

9627E-01 

0 

1082 

44 

00 

-0 

9516 

-1  050 

0 

9857E-01 

0 

7024 

45 

00 

-0 

9775 

-1.  077 

0 

99S1E-01 

0 

2969 

46 

00 

-0 

9927 

—1  094 

0 

9997E-01 

0 

2915 

47 

00 

-0 

9999 

-1.  099 

0 

9S91E-01 

0 

2867 

43 

00 

-0 

9962 

-1  092 

0 

964SE-01 

0 

2312 

49 

00 

-0 

9325 

-1.  075 

0 

9255E-01 

0 

2765 

50 

00 

-0 

9539 

-1  046 

0 

8697E-01 

0 

2718 

51 

00 

-0. 

9258 

-1  005 

0. 

7962E-01 

0 

2672 

52 

00 

-0 

3325 

-0  9523 

0 

7029E-01 

0 

2629 

c T* 

-_i 

00 

-0 

8222 

-0  8914 

0 

5915E-01 

6i 

2587 

54 

00 

-0 

7723 

-0  8136 

0 

45:2£iE_01 

0 

2546 

55 

00 

-0. 

7055 

-0.  7253 

0 

2025E-01 

6i 

25£i7 

56 

00 

-0 

6212 

-0.  6427 

0 

1241E- 01 

6i 

2463 

57 

00 

-0 

5507 

-0  5429 

-0 

7760E-02 

£i 

2421 

58 

00 

-0 

4646 

-0  4247 

-0 

2022E-01 

0 

2 2' 9 5 

5? 

00 

-0 

-?7“  q 

-0.  2136 

-0 

552SE-01 

i? 

2260 

60 

00 

-0 

2794 

-0  1968 

-0 

S265E-01 

0 

61 

00 

-0 

1822 

-0  6975E-01 

-0 

1124 

0 

2292 

Fig.  3 - Program  ADAPT  Sample  Output 
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0.9857 

0.05714 

-0.0857 

0.05714 

-0.01429 

1.488 

-1.619 

-0.5714 

1.048 

-0.  1452 

0.6429 

-1.071 

-0.1429 

0.9286 

-0.3571  , 

0.08333 

-0.1667 

0.0 

0.1667 

-0.08333  i 

Weight  Matrix 


— 

0.9857 

1.488 

0.6429 

0.08333 

1.488 

6.379 

3.869 

0.5972 

1 

| 

0.6429 

3.869 

2.571 

0.4167  ! 

i 

0.08333 

0.5972 

0.4167 

0.06944 

i 

Covariance 

Matrix 

Table  1 Optimal  Weight  Matrix  and  Covariance  Matrix  For 

a Five  Point  Nonrecurslve  Filter  Using  a Third 
Order  Polynomial  Model.  The  Filter  Estimates 
the  Data  at  the  End  of  the  Window 
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-0.08571 

0.  1429 

0.4857 

O.  1470 

-0.085/1 

-0.08333 

0.6667 

0.0 

-0.666/ 

0. OH  111 

0.1429 

-0.07143 

-0.1429 

-0.07143 

0.1429 

0.08333 

-0.1667 

0.0 

0.1667 

-0.08333 

Weight 

Matrix 

0.4857 

0.0 

-0.1429 

0.0 

0.0 

0.9828 

0.0 

-0.2361 

-0.1429 

0.0 

0.07143 

0.0 

0.0 

-0.2  16 1 

0.0 

0.06944 

Covarianre  Matrix 


Table  2 Optimal  Weight  Matrix  and  Covariance  Matrix  For 
a Five  Point  Nonrecursive  Filter  Using  a Third 
Order  Polynomial  Model.  The  Filter  Estimates 
the  Data  in  the  Center  of  the  Window 


i 
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INTRODUCTION 

This  report  contains  the  documentation  for  the  LPASS  program. 

It  consists  of  the  design  procedure  used,  a description  of  the 
program,  and  design  examples  using  the  program. 

The  purpose  of  the  LPASS  program  is  the  design  of  a maximally 
flat  Butterworth  or  an  equiripple  Chebychev  lowpass  digital  filter. 
Starting  with  an  analog  filter,  the  bilinear  Z transform  is  used  to 
find  an  equivalent  digital  filter.  The  user  enters  the  following 
parameters:  the  number  of  second  order  sections,  the  type  of  filter, 

the  sampling  interval,  the  -3db  cutoff  frequency,  the  starting  frequency 
and  the  frequency  increment.  If  a Chebychev  filter  is  being  designed, 
the  ripple  must  also  be  entered. 

The  program  calculates  the  digital  filter  coefficients  for  up  to 
three  second  order  sections  in  cascade.  The  program  is  designed  to 
calculate  up  to  a sixth  order  filter,  thus  the  filter  order  is  two 
times  the  number  of  cascaded  second  order  sections.  The  filter 
magnitude  response  is  generated  over  the  frequency  interval  specified 
by  the  input. 

The  LPASS  program,  written  in  Fortran  IV,  is  supplied  as  a card 
deck  with  this  report.  The  program  is  in  the  form  of  a subroutine 
and  can  be  used  as  is  by  a call  statement  from  the  main  program. 

Data  may  be  input  via  cards  with  output  available  through  a line 
printer.  The  input/output  devices  may  be  altered  as  explained  in 
this  report.  Graphics  routines  may  easily  be  appended  to  the  program. 
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I . Design  Procedure 

A.  Preliminary  Discussion 

The  transfer  function  of  a second  order  digital  filter  in  the 
Z domain  is  given  by 


H(Z) 


K1(AQZi+A1Z+A2) 


Z“+81Z+B2 


(1) 


where  the  A's  and  B's  are  the  coefficients  of  the  numerator  and 
denominator  respectively.  One  common  method  of  designing  a digital 
filter  is  to  start  with  an  analog  transfer  function  H(S)  and 
transform  it  to  the  digital  transfer  function  H(Z).  This 
program  will  calculate  the  scale  factor  Kj  and  the  coefficients 
Aq,  A^,  A2>  B^,  and  B2<  The  transformation  used  is  the  extended 
bilinear  Z transform  defined  as 


2 (Zzij 

t vz+r  ’ 


(2) 


where  T is  the  sampling  interval.  When  this  transform  is  employed 
the  desired  frequencies  must  first  be  prewarped  to  make  them  com- 
patible with  the  digital  filter.  The  prewarped  cutoff  frequency  is 
given  by 


WDC 


u T 


(3) 


This  prewarping  is  done  by  the  program. 

B.  Butterworth  Low-Pass  Filter 

We  start  with  a normalized  second  order  low-pass  filter  in 
the  S plane. 


1 


H(S)  = 


2 

S + 2Sc.os0  + 1 


O) 
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where  the  angle  0 is  in  degrees  (in  the  program).  0 may  be  found 
from  the  Butterworth  circle  and  the  relationship 

g = e±jir(2m-l)/2n  (5) 

where  n is  the  order  of  the  filter  and  m = 1,  2,  3,  ...  , n. 

This  relationship  is  determined  by  the  following  procedure.  By 
definition,  a filter  is  n*"'1  order  Butterworth  lov-pass  if  its  gain 
characteristic  is 

_ 2 


|Hn(ju)) 


(6) 


1 + (— ) 
to 

c 


2n 


where  a is  the  gain,  loc  is  the  desired  cutoff  frequency  and  n is 

i 1 2 

the  order  of  the  filter.  Note  that  j ( j co) [ goes  to  zero  as  to  goes 
to  infinity,  indicating  the  filter  does  attenuate  the  higher  frequencies. 
To  determine  its  efficiency  as  a low-pass  filter  we  calculate 


— |H  (j<o)  | = -•  — 
d<o'  n 1 to 

c 


2n- 

(— ) 

0) 

c 

■1 

r 2n" 

wr> 

L c J 

3/2 

(7) 


Thus 


_d 

dio 


H (jio) 
ri 


(8) 


to=0 


for  all  n and  hence  the  gain  characteristic  stays  flat  for  co  close 
to  0.  Also 


-t-|h  (» 
dto 1 n J 


an 


<o=io  2io  /jf 

c C 


(9) 


and  hence,  the  decline  rate  or  "roll-off"  of  the  gain  characteristic 


J 


at  0)  = u becomes  sharper  as  n increases.  In  other  words,  the 
approximation  to  the  ideal  low-pass  filter  improves  for  larger  n. 

The  order  n is  chosen  according  to  desired  specifications.  The 
references  have  equations,  curves,  and  tables  that  select  n,  given 
the  specifications.  For  example,  page  227  of  Rabiner  and  Gold  gives 
an  equation  for  calculating  n when  the  transition  band  is  specified. 

In  the  design,  the  poles  for  the  full  frequency  response,  H(S), 
of  the  n*"*1  order  Bucterworth  filter  must  be  determined.  The  pro- 
cedure is  as  follows: 


,2 


I Hn  (jco)  j = Hn(ju))Hn(jw)  = Hn(jw)Hn(jw)  = Hn(ju))Hn(-jw) 


= [H(S)H(-S) ] 


S=ju> 


2n 


1 + (-^) 


2n 


CO  CO  = T 

c 3 


S 1 + (VH 


1 + - 


2 

•-w 

c 


1 + 


for  n even 


c J 


1 f 
1 - |~! 

K j 


, for  n odd 


(10) 


Setting  the  denominators  equal  to  zero, 


_S 

CO 

c 


(il) 


l/2n 


(ID 


Thus,  the  pole  locations  are  the  2n  roots  of  ±1,  depending  on  whether 
the  order  is  odd  or  even.  These  roots  are  located  on  a circle  with 
radius  centered  at  the  origin  of  the  S plane  and  have  symmetry 
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with  respect  to  both  real  and  imaginary  axes.  For  n odd,  a pair  of 
roots  are  on  the  real  axis  and  the  rest  are  separated  by  u/n  radians. 
For  n even,  a pair  of  roots  are  located  n/2n  radians  from  the  real 
axis  and  the  rest  are  again  separated  by  ir/n  radians.  No  roots  are 
on  the  imaginary  axis  for  either  even  or  odd  n. 


Let  p^,  • • • » 
locations,  if  p^,  . . 


?2n  be  the  roots.  From  the  symmetry  of  the  pole 
. , pn  are  the  roots  lying  in  the  right-half 


plane,  the  left-half  plane  roots  are  -p^,  • 


-p  . The 
rn 


magnitude-squared  function  can  then  be  written  as 
H (S)H  (-S)  = 


2 , , . n „n 
a (-1)  a)  2 
c 


(s+Pl)...(s+pn)(s-Pl)...(s-pn) 


To  bo  stable,  H (S)  must  have  all  its  poles  In  the  left-half  piano,  thus 


aw 


Hn(S)  (S+Pl)...(S+pn) 


(13) 


The  program  is  written  with  unity  gain  at  DC,(co=0),  therefore  a = 1. 

In  order  to  locate  the  poles  as  specified  above,  consider  the 
following  set  of  equations. 


. ±j  tt  (2m-l) 

1 = -e  , m = 1 , 2 , . . . , n;  for  n even 


-1  = -e 


±j  2-rrk 


(14) 


, k = 0,  1,  . . . , n;  for  n odd 


Substituting  equations  (14)  into  equation  (11)  yields 
r_Si  _ ±j*(2m-l)/2n  _ , 

co  -‘+m  ~ -e  , m - 1,  2,  . . . , n;  for  n even 

c 

(15) 

r S-,  +jirk/n  , 

co  +k  _ -e  , k = 0,  1,  . . . , n;  for  n odd 

c 

Equations  (15)  will  give  the  pole  locations  as  described  above. 
Consider  the  form  of  equations  (15) 

+ 10 

S = -co  e =u)  f-cos0  ± jsinGl. 
c c 1 


(16) 
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From  this  relationship,  it  can  be  seen  that  the  magnitude  for  each  pole 
is  regardless  of  the  angle,  and  thus  all  the  poles  lie  on  a circle 

with  radius  w . 

c 

As  an  example  consider  a second  order  filter,  n = 2. 
r St  _ ±jir(2m-l)/ 4 m = 1,  2 

[- — J+™  -e 

to  ±m 
c 

S , , = w / ±45° 

±1  c 

Si2  = u / ±133° 


0 = 45° 


The  relationship  of  these  roots  about  the  circle  of  radius 
is  illustrated  in  Figure  1.  The  angle  0 is  always  measured  from 
the  negative  real  axis. 

In  the  program,  only  the  angle (s)  less  than  90°  are  considered 
so  that  poles  lie  in  the  left-half  plane  since  poles  in  the  left-half 
plane  are  stable.  Putting  0 = 45°  into  equation  (4)  yields  poles 
at  -0.707  ± J0.707.  These  locations  are  in  the  left-half  plane. 

In  the  program,  only  even  order  filters  are  considered. 

Below  are  the  values  of  0 for  1,  2,  and  3 second  order  sections 
in  cascade. 


Cascaded 

Sections 


Filter 

Order 


Angle 


N 


n 0 


1 

2 

3 


2 45° 

4 22.5°,  67.5° 

6 75°,  45°,  15° 


These  calculated  angles  are  incorporated  in  the  program  in  the  order 


given  above. 


For  N second  order  sections  there  are  N 0's.  Only  one  specific 
0 is  used  per  stage,  because  each  stage  has  only  one  set  of  pole 
locations. 

The  following  is  the  procedure  to  derive  the  magnitude  of  the  ith 
stage,  where  i varies  from  1 to  N. 

Given  the  normalized  second  order  low-pass  transfer  function 
equation  (4),  we  employ  the  low-pass  to  low-pass  transformation  for 
an  arbitrary  cutoff  frequency  u>c  given  by 


S 


(17) 


til 

For  the  i stage,  equation  (4)  becomes 


2 

H.(S)  = -r J (18) 

S 42Sw  cos©  ,+uj 
c i c 


The  extended  bilinear  Z transform,  equation  (2),  is  used  to  get  to 


the  digital  domain.  Employing  equation  (2)  on  equation  (18)  and  sub- 


stituting WDC  for  yields 


Hi  (Z) 


MDC2(Z2+2Z+1) 

-|(Z2-2Z+1)+|(Z2-1)WDC  cos©i+WDC2(Z2+2Z+l) 


(19) 


Putting  the  denominator  of  equation  (19)  in  monic  form  yields  the 
transfer  function  for  the  i*"^1  stage  of  the  filter 


H^Z)  = 


k11(aqz>a1z+a2) 

Z2+Bli2+B2i 


(20) 


Equation  (20)  is  the  same  as  equation  (1)  with  the  exception  of  the 
subscripts.  In  equation  (20) 
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Gi  = -|  + | WOC  cos0i  + WDC2 


(21) 


11 


11 


21 


WDC 


2 WDC2  - 


-j  - ^ WDCcose  +WDC2 
T2  1 1 


ST 

Letting  Z = e and  S = joa  and  taking  the  magnitude  of  H^Cjw)  we 
have 

J( A cos  (2uT)+A  cos (a>T)+A  ) 2+  (A  sin (2wT)+A  sin (wT)  ) 2 

|h  a«)i  = K — ■ ■ - S r 

J (cos(2u)T)+B11co8(a)T)+B2i)  +(Sin(2a)T)+B1;.sin(uT)) 


(22) 


This  magnitude  function  is  the  same  for  both  the  Butterworth  and  the 
Chebychev  filters  where  i varies  from  1 to  N. 

C.  Chebychev  Low-Pass  Filter 

The  advantage  of  the  Chebychev  low-pass  filter  over  the 
Butterworth  low-pass  filter  is  that  the  transition  band  of  the 
response  at  frequencies  greater  than  is  sharper  for  the  Chebychev 

low-pass  filter.  This  is  achieved  by  specifying  a small  percentage 
of  ripple  in  the  low-pass  region.  The  amplitude  of  the  ripple  is 
specified  by  the  quantity  6 (labeled  RIP  in  the  program).  Figures  6, 
7 and  8 illustrate  the  rippling  for  second,  fourth,  and  sixth  order  fi- 
ters,  respectively.  The  poles  of  the  filter  are  found  on  an  ellipse 
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described  by  two  Butterworth  circles  of  radii  A and  B with  A<B. 
The  location  of  the  poles  on  the  ellipse  is  a function  of  the  ripple, 
6,  and  is  given  by  the  following  equation: 

-1.-1/2N, 


B,  A = j((Je~2+l +e_1)1/2N  ± (JeF-l+e  x)  x' x") 


(23) 


where 


1 

<i-«F 


i 


1/2 


(24) 


B is  given  for  the  plus  sign  and  A for  the  minus  sign.  The  Chebychev 
ellipse  then  has  major  axis  B and  minor  axis  A.  The  location  of  the 
S plane  poles  on  the  ellipse  is  given  by 


Real  Part  = A cos0 


Imaginary  Part  = B sin0 


(25) 


The  0’s  are  the  same  as  given  for  the  corresponding  order  Butterworth 
filter.  An  example  of  Chebychev  pole  locations  is  illustrated  in 
Figure  2.  For  A ■ 1/2  and  B = 1 in  a fourth  order  filter, 

0 = 22.5°  and  67.5°.  The  Chebychev  pole  locations  are  determined 
from  equations  (23),  (24),  and  (25). 

The  analog  second  order  Chebychev  low-pass  filter  is 


H(S) 


K2a 


S +KgS+K2 


(26) 


where 


a = 


1/N 


(27) 


e is  calculated  from  equation  (24)  and  N is  the  number  of  second 
order  sections.  Kg  and  K2  are  calculated  by 
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Kg  = 2Acos6 


2 2 2 2 
K2  = A cos  0 + B sin  0 . 


The  substitution  of  the  low-pass  to  low-pass  transformation  ror  some 
cutoff  frequency  oj^,  equation  (17),  into  equation  (26)  yields 


H (S  / = 


K-aoj 

l c 

2 2 
S + SK^O)  + u K„ 
8 c c 2 


Using  the  extended  bilinear  Z transform,  equation  (2),  and  subsr<n.Mno 

WDC  for  uj  we  have  for  any  section 
c 

K aWDC2(Z2+2Z+l) 

u,7x 2 . (31) 

H(Z)  ~/o  9 9 2 2 k 

— |(Z  -2Z+1)  + |KgWDC(ZZ-l)  + WD(TK2(Z  +2Z+1) 


Collecting  terms  yields  the  following  for  the  i section 


H±(Z) 


Kli(A0Z  +A1Z+A2) 


Z +BliZ+B2L 


where 


A0  = a2  1 


Al=  2 


Gi  = ~1  + I WDC- Kg  + WDC2K2 
aK2WDC2 


2WDC2K  - -| 
2 Tz 
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~~  - | WDC-Kg  + WDC2K2 

B2i  ” 

i varies  from  1 to  N.  These  coefficients  are  used  to  find  |H^(jw)| 
given  in  (22). 

For  applications  where  a sharper  roll-off  is  required  the  Chebychev 
filters  are  used.  The  roll-off  increases  with  n for  any  fixed  e. 

For  fixed  n,  the  roll-off  decreases  as  e decreases.  For  small 
e the  ripple  width,  6,  is  small,  see  equation  (23),  but  so  is  the 
roll-off.  For  larger  e the  roll-off  improves  but  the  ripple  width 
increases.  In  the  first  case  the  filter  will  be  good  at  L/C  ana 
low  frequencies,  unsatisfactory  at  high  frequencies.  The  converse  is 
true  in  the  second  case. 

The  above  observations  suggest  the  procedure  to  be  used  in 
selecting  a Chebychev  filter  to  match  a set  of  specifications.  The 
permissible  ripple  width  specifies  e.  With  e fixed,  select  n 
to  attain  the  required  roll-off. 


r 
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II.  Using  the  Program 

The  first  data  card  read  into  the  program  contains  the  number  of 

second  order  sections  to  be  cascaded,  N,  and  the  type  of  filter 

desired,  KN.  N is  equal  to  1,  2,  or  3,  which  corresponds  to 

the  2nC*,  4t*1,  or  6^  order  filter  respectively.  KN  = 1 yields  a 

Butterworth  filter,  while  KN  = 2 yields  a Chebychev  filter.  The 

format  on  the  N,  KN  card  is  212.  The  second  data  card  read  in 

is  the  sampling  interval  T in  F10.6  format.  When  choosing  T,  1/T 

should  be  approximately  equal  to  ten  times  the  cutoff  frequency,  w . 

The  third  data  card  contains  the  value  of  oj  in  F10.4  format. 

c 

For  the  Butterworth  low-pass  filter,  is  the  -3db  cutoff  frequency. 

2.1/2 

For  the  Chebychev  filter  the  magnitude  of  the  response  is  i/(l+c  ) 

= 1-6  at  a)  = a)  . to  is  in  radians.  6 is  the  ripple  factor. 

If  the  desired  filter  is  Chebychev,  i.e.,  KN  = 2,  the  next 

data  card  is  the  ripple  factor  (RIP)  in  F5.3  format.  The  filter 

response  for  all  even  order  Chebychev  low-pass  filters  passes  through 
2 1/2 

l/(l+e  ) =1-6  for  u = 0 and  w . For  odd  order  filters,  the 

c 

2 1/2 

magnitude  is  1 for  w = 0 and  l/(l+e  ) =1-6  for  u = 

Tills  program  produces  only  even  order  filters.  If  the  desired 
filter  is  Butterworth,  i.e.,  KN  = 1,  this  data  card  is  omitted  from 
the  data  deck. 

The  final  data  card  is  the  starting  frequency  (FREQ1)  and  the 

frequency  increments  (DELT)  in  radians.  The  format  of  the  FREQ1, 

DELT  card  is  2F10.4.  Determine  DELT  by  the  following: 

D T _ final  frequency  - starting  frequency 

1024 


This  is  necessary  because  there  are  1024  frequency  data  points  calculated 
in  the  program.  Choose  FREQ1  and  DELT  to  insure  that  calculated  values 
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will  Include  the  data  of  interest.  For  maximum  efficiency  of  the 
program,  DELT  should  be  a multiple  of  2 so  no  decimal  to  binary 
conversion  errors  are  incurred. 

The  digital  filter  coefficients  are  computed  and  printed  out  for 
each  second  order  section.  The  full  filter  magnitude  response,  as  well 
as  each  section  magnitude  response,  is  printed  for  each  of  the  frequency 
increments  specified.  When  N = 1,  the  section  magnitude  response  is 
the  full  filter  magnitude  response  and  is  only  printed  once. 

The  program  may  be  easily  modified  to  incorporate  a graphics 
display  of  the  magnitude  response.  There  is  a comment  card  in  the 
LPASS  program  indicating  where  the  graphics  subroutine  call  card  should 
be  inserted. 

The  program  is  written  with  input  obtained  via  device  4 and  output 
written  to  device  6.  These  numbers  should  be  assigned  to  the  appro- 
priate devices  prior  to  running  the  program. 

The  program  was  developed  on  the  PDP-11/20  with  a DOS/BATCH 
operating  system.  Trial  runs  frequently  used  a TTY  terminal  as 
well  as  a card  reader  for  input  (device  4);  and  a TTY  terminal  as 
well  as  a line  printer  for  output  (device  6) . 

Double  precision  arithmetic  is  employed.  To  decrease  required 
memory  storage,  only  the  frequency  interval  values  and  the  full 
magnitude  response  are  saved.  The  section  magnitude  responses  are 
printed  out,  but  are  not  stored.  The  program  will  produce  approxi- 
mately 21  pages  of  output. 

Shown  below  are  sample  deck  set-ups  for  the  Chebychev  and 
Butterworth  low-pass  filters. 
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i 


Data 

Card 

Format 

Example 

1 

212 

0302  (3  sections  Chebychev  low-pass) 

2 

F10.6 

0.001  (T  = 0.001) 

3 

F10.4 

100  (to  = 100  radians) 
c 

4 

F5.3 

0.10  (Ripple  amplitude  = 0.10) 

5 

2F10.4 

70  0.06  (Start  at  w = 70.  Steps  oi 
0.06  radians.  Will  finish  just 
past  w = 131  radians.) 

1 

212 

0201  (2  sections,  4C^  order,  Butterw. 

2 

F10.6 

0.005  (T  = 0.005) 

3 

F10.4 

20  (w  = 20  radians) 
c 

4 

2F10.4 

0 0.04  (Start  i*  j 0 

past  w = 40  radians  in  steps  of 

0.04  radians.) 


The  following  pages  contain  annotated  examples  of  output  data. 
This  is  an  example  of  the  output  for  a 41*1  ofder  Buttervorth 
low-pass  filter  with  T » 0.005  and  = 20  radians.  The  starting 
frequency  is  0 radians  and  the  frequency  increment  is  0.04  radian. 
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I is  the  stage.  I varies  from  1 to  N. 

WDC  is  the  prewarped  cutoff  frequency. 

WC  is  the  cutoff  frequency. 

T is  the  sampling  interval. 

Aq , A^,  and  are  the  low-pass  filter  numerator  coefficients, 
and  B2  are  the  low-pass  filter  denominator  coefficients, 
is  the  gain  factor. 

W is  the  frequency. 

H is  the  overall  magnitude  of  the  digital  transfer  function. 

s t 

HI  is  the  magnitude  of  the  digital  transfer  function  (1  stage). 

H2  is  the  magnitude  of  the  digital  transfer  function  (2n<^  stage). 

H = H1*H2 . 

See  Figure  4. 

This  1 s an  example  of  the  output  for  a 6t'1  order  Chebychev 
low-pass  filter  (three  second  order  stages  cascaded)  with  T = 0.005 
and  wc  = 20  radians.  The  starting  frequency  is  0 and  the  frequency 
increment  is  0.04  radian.  The  ripple  is  equal  to  0.100. 


WDC  - 20.01668  WC  = 20.00000  T = 

A = 0.24783947  B = 1.03025433 

A = 0.24783947  B = 1.03025453 

A = 0.24783947  B = 1.03025453 

FOR  I = 1 AQ  = 0. 10000000E+01 
Ax  = 0.10000000E+01 
B1  = -0. 19774006E+01 

WDC2  = 0.40066761E+03  G(I)  - 0. 

FOR  I = 2 AQ  = 0.10000000E+01 
A2  = 0. 10000000E+01 
B1  = -0. 19600541E+01 


0. 50000E-02 

= 0.12829114  K = 0.99443709 
1C  = 0.35049793  KZ  = 0.56142438 
Kg  = 0.47878908  = 0.12841170 

A1  = 0.20000000E+01 
K1  = 0. 23830688E-02 
B2  = 0. 98727357E+00 

16142562E+06  A = 0. 96548939E+00 

Ax  = 0.20000000E+01 

^ = 0. 13321469E-02 

E = 0.96557320E+00 
2 
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WDC2  = 0. 40066761E+03 
FOR  I = 3 


G(I)  = 0.16303127E+06  A = 0. 96548939E+00 


Aq  = 0. 10000000E+01 
A2  = 0. 10000000E+01 
B1  = -0. 19519613E+01 


Ax  = 0. 20000000E+01 
K1  = 0. 30310788E-03 
B2  = 0. 95321709E+00 


WDC2  = 0. 40066761E+03  G(I)  - 0. 16388496E+06  A = 0. 96548939E+00 


w 

H 

III 

H2 

Hi 

0.0000 

0. 89998E+00 

0. 96549E+00 

0.96549E+00 

0. 96547E+00 

0.0400 

0. 89999E+00 

0. 96549E+00 

0. 96549E+00 

0.96548E+00 

0.0800 

0. 90001E+00 

0. 96550E+00 

0. 96551E+00 

0. 96547E+00 

0. 1200 

0. 90009E+00 

0. 96552E+00 

0.96554E+00 

0. 96550E+00 

0.1600 

0. 90016E+00 

0.96555E+00 

0. 96558E+00 

0. 96551E+00 

0.2000 

0. 90028E+00 

0. 96558E+00 

0.96564E+00 

0. 96555E+00 

0.2400 

0. 90042E+00 

0. 96563E+00 

0. 96571E+00 

0.  9605^1.-1-0^ 

• 

. 

• 

• 

• 

• 

• 

• 

. 

WDC  is  the  prewarped  cutoff  frequency. 

WC  is  the  cutoff  frequency. 

T is  the  sampling  interval. 

b,  a - i (<J7W1)1/2N!<1/rW-V1/2N> 

I is  the  i1"^1  stage,  I varies  from  1 to  N. 

Kg  = 2Acos0. 

2 2 2 2 
K2  = A cos  6+B  sin  0. 

Aq,  Aj,  and  A2  are  the  low-pass  filter  numerator  coef f iciei.ts. 
B^  and  B2  are  the  low-pass  filter  denominator  coefficients. 

K^  is  the  gain  factor. 

WDC2  = (WDC)2. 

G(I)  = + |wDC-Kg  + (WDC)2K2. 


The  A following  C(I)  is  a 


.1 


1 11/2N 

77, 


W is  the  frequency. 

H is  the  overall  magnitude  of  the  digital  transfer  function. 
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H1  Is  the  magnitude  of  the  digital 
H2  is  the  magnitude  of  the  digital 
H3  is  the  magnitude  of  the  digital 
H = H1*H2*H3. 


transfer 

function 

dSt 

stage) 

transfer 

function 

(2nd 

stage) 

transfer 

f unction 

(3rd 

stage) 

See  Figure  8. 
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INTRODUCTION 

This  report  contains  the  documentation  for  the  BPASS  program.  It 
consists  of  the  design  procedure  used,  a description  of  the  program, 
and  design  examples  using  the  program. 

The  purpose  of  the  BPASS  program  is  the  design  of  either  a 
maximally  flat  Butterworth  or  a Chebychev  filter  with  equal  ripple  in 
the  pass  band.  For  each  type  of  filter  there  is  a choice  of  band-pass 
or  band-stop  filters.  Starting  with  an  analog  filter,  the  bilinear 
Z transform  is  used  to  design  an  equivalent  digital  filter.  The  user 
enters  the  low-pass  filter  order,  the  type  of  filter  desired,  the 
sampling  interval,  the  upper  and  lower  cutoff  frequencies,  the 
starting  frequency  and  frequency  increment,  and  if  a Chebychev  filter 
is  being  designed,  the  ripple.  The  low-pass  filter  sections  are 
transformed  to  second  order  band-pass  or  band-stop  sections.  Then  the 
program  generates  the  digital  filter  coefficients  for  up  to  six 
second  order  sections  in  cascade  or  up  to  a 12th  order  filter.  The 
design  is  carried  out  in  the  frequency  domain.  The  program  calculates 
the  transfer  function  coefficients  for  each  second  order  section,  the 
magnitude  function  for  each  section,  and  the  final  cascaded  filter 
magnitude  response  over  the  frequency  interval  specified  by  the  input. 

The  BPASS  program,  written  in  Fortran  IV  is  supplied  as  a card 
deck  with  this  report.  The  program  is  in  the  form  of  a subroutine  and 
can  be  used  as  is  by  a call  statement  from  the  main  program.  Data 
may  be  input  via  cards  with  output  available  through  a line  printer. 
The  input /output  devices  may  be  altered  as  explained  in  this  report. 


Graphic  routines  may  easily  be  appended  to  the  program. 


r 

* 
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I . Design  Procedure 

A.  Preliminary  Discussion 

One  common  method  of  designing  a digital  filter  is  to  start 
with  an  analog  transfer  function  H(S)  and  transform  it  to  the 
digital  transfer  function  H(Z). 

The  transfer  function  of  a second  order  digital  filter  in  the 
Z domain  is  given  by 

K. (AnZ2  + A Z + A ) 

H(Z)  = i — (1) 

Z + BjZ  + B2 

where  the  A's  and  B's  are  the  coefficients  of  the  numerator  and 
denominator  respectively.  This  program  will  calculate  the  scale 
factor  and  the  coefficients  Aq,  A^,  A2>  B^,  and  The 

transformation  used  is  the  extended  bilinear  Z transform 


S 


2.Z  - 1 
T 'Z  + 1 


) 


(2) 


where  T is  the  sampling  interval.  When  the  extended  bilinear  Z 
transform  is  employed,  the  desired  frequencies  must  first  be  pre- 
warped to  make  them  compatible  with  the  digital  filter.  In  the 
band-pass  and  band-stop  filters,  the  upper  and  lower  cutoff  frequencies 
and  the  center  frequency  of  the  filter  are  of  interest.  Calling  the 
upper  and  lower  frequencies  and  respectively,  the  pre- 

warped upper  (WDU) , lower  (WDL) , and  center  (WDM)  frequencies  and  the 
bandwidth  between  WDU  and  WDL  are  found  by 
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2 W T 
WDU  = - tan(-y-) 

2 

WDL  = - tan(— ) 


„ no  u)  T 

WDM  = ~ tan  — ~~ 


WB  = WDU  - WDL 


0)^  and  <o^  are  specified  by  the  designer  and  the  prewarping  is  done 
by  the  program. 

In  the  design  procedure  for  all  band-pass  and  band-stop  filters  of 
order  n*.  (n'  even),  the  program  begins  by  first  finding  the  poles 

for  the  corresponding  n'/2  order  low-pass  filter.  The  low-pass 
filter  is  then  transformed  into  a band-pass  or  band-stop  filter  of 
order  n’,  (n'  = 2n). 

B.  Butterworth  Band-Pass  Filter 

We  start  with  a normalized  second  order  low-pass  Butterworth 
filter  transfer  function  in  the  S plane 


S + 2Scos0  + 1 

where  the  angle  0 is  in  degrees  (in  the  program)  and  may  be  found 
from  the  Butterworth  circle  and  the  relationship 


S » e 


±jtr(2m  - l)/2n 


where  n is  the  order  of  the  low-pass  filter  and  m ■ 1,  2,...,  n. 
This  relationship  is  determined  by  the  following  procedure.  By 


definition,  a filter  is  nth  order  Butterworth  low-pass  if  its  gain 


characteristic  is 
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|Hn(» 


c 


2n 


(6) 


where  a is  the  DC  gain,  w is  the  desired  cutoff  frequency  and  n 

c 

is  the  order  of  the  low-pass  filter. 

In  the  design,  the  poles  of  H(S)  must  be  found.  The  procedure 
is  as  follows: 


|Hn(jU)r  = Hn(jW)Hn(ju)  = Hn(ju)Hn(jU)  = Hn(ju)Hn(-ju) 


= [H(S)H(-S) ] 


S = ju 


1 + (*-) 
0) 


2n 


c JO)  = 


2n 


s 


k" 

n 

1 + 

2 

0) 

c 

, for  n even 


1 + 


.2  i 


c J 


1 - 


rslt 

-2 

U) 

L c J 


, for  n odd 


Setting  the  denominators  equal  to  zero, 


S_  . (±i)l/2n 
0) 

c 


(7) 


(8) 


Thus,  the  pole  locations  are  the  2n  roots  of  +1,  depending  on 
whether  the  low-pass  filter  order  is  odd  or  even.  These  roots  are 
located  on  a circle  with  radius  u)c  centered  at  the  origin  of  the  S 
plane  and  have  symmetry  with  respect  to  both  real  and  imaginary  axes. 
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For  n odd,  a pair  of  roots  are  on  the  real  axis  and  the  rest  are 
separated  by  n/n  radians.  For  n even,  a pair  of  roots  are  located 
ir/2n  radians  from  the  real  axis  and  the  rest  are  again  separated  by 
ir/n  radians.  No  roots  are  on  the  imaginary  axis,  for  either  even  or 
odd  n. 


Let  p ,...,p2n  be  the  roots.  From  the  symmetry  of  the  pole 
locations,  if  p ,...,pn  are  the  roots  lying  in  the  right -half  plane, 
the  left-half  plane  roots  are  -p^,...,-pn.  The  magnitude-squared 
function  can  then  be  written  as 


2.  .n  2n 
a (-1)  wc 

Hn(S)Hn(_S)  = (S  + Pl)...(S  + pn)(S  - P;L)...(S  - pn)  * 


(9) 


To  be  stable,  H (S)  must  have  all  its  poles  in  the  left-hand  plane, 
thus 


Hn(S) 


aoj 


(S  + P . ) • . . (S  + p ) 
J.  n 


(10) 


The  program  is  written  with  unity  gain  at  DC,  (w  - 0) , therefore 
a * 1. 

In  order  to  locate  the  poles  as  specified  above,  consider  the 
following  set  of  equations. 

1 = ^ , m = 1,  2,...,  n;  tor  n even 

-1  = -e1^2^  , k = 0,  1,...,  n;  for  n odd 


(ID 


Substituting  equations  (11)  into  equations  (8)  yields 


lu_ )+in 
c 


j In (2m  - l)/2n  , , c 

-e  J , m-l,2,..,n;  for  n even 


(12) 


ri_i  = 
lu  J±k 
c 


, k ■ 0,  1,...,  n;  for  n odd 
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Equations  (12)  will  give  the  pole  locations  as  described  above. 
Consider  the  form  of  equations  (12) 


S = -w  w [-cos0  + jsine]  . (13) 

c c 

From  this  relationship,  it  can  been  seen  that  the  magnitude  for  each 

pole  is  regardless  of  the  angle,  and  thus  all  the  poles  lie  on  a 

circle  with  radius  w . 

c 

As  an  example,  consider  a second  order  filter,  n = 2. 


r S_i  _ ±jir(2m  - l)/4 

J±m  “ 

c 


m = 1 , 2 


S 


±1 


= u>  /+ 45° 
c 


S.,  = w /+135° 

±1  c 

0 - 45° 

The  relationship  of  these  roots  about  the  circle  of  radius  is 
illustrated  in  Figure  1.  The  angle  0 is  always  measured  from  the 
negative  real  axis. 

In  the  program,  only  the  angle (s)  less  than  90°  are  considered 
so  that  the  poles  lie  in  the  left-half  plane  because  poles  in  the 
left-half  plane  are  stable.  Putting  0 ■ 45°  into  equation  (4) 
yields  poles  at  -0.707  +J0.707.  These  locations  are  in  the  left-half 
plane.  From  equations  (12),  for  low-pass  filter  orders  n ■ 1,  2 


6,  the  values  of  0 are  given  below. 
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Second 

Band-Pass 

.ow-Pass 

Order 

Band-Stop 

Filter 

Cascaded 

Filter 

Order 

Angle 

Sections 

Order 

n 

0 

N 

n' 

1 

0° 

1 

2 

2 

45° 

2 

4 

3 

60°,  0° 

3 

6 

4 

22.5°,  67.5° 

4 

8 

5 

72°,  36°,  0° 

5 

10 

6 

75°,  45°,  15° 

6 

12 

n is  the  order  of  the  low-pass  filter  and  is  used  to  determine  pole 
locations,  n is  also  the  number  of  second  order  band-pass  or  band- 
stop  sections  which  results  from  the  transformation  of  the  low-pass 
filter  sections  and  which  will  be  cascaded  to  form  the  band-pass  or 
band-stop  filters  of  order  n'.  The  transformation  is  explained  below. 
The  calculated  angles  are  incorporated  in  the  program  in  the  order 
given  above. 

Given  the  normalized  second  order  low-pass  transfer  function 
equation  (4),  we  transform  this  low-pass  into  a band-pass  transfer 
function  for  some  bandwidth  WB,  and  center  frequency  WDM  by 
us’n'  the  transform 


2 2 
S + WDM 

S -*■  

SWB 

Equation  (4)  then  transforms  to  a 4th  order  transfer  function 
H(S)  = 


(14) 


2 2 
S WB 


S4  + S3  2WBcos0  + S2 (2 WDM2  + WB2)  + S2WB  WDM2  cos6  + WDM4 


(15) 


Using  the  root  finding  subroutine  "POLRT"  from  the  IBM  Scientific 
Subroutine  Package  (SSP),  the  roots  of  the  denominator  of  equation  (15) 


p 
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are  found.  (Note:  POLRT  has  been  attached  to  BPASS  as  a double 
precision  subroutine  and  is  included  in  the  card  deck).  The  roots 
found  will  be  coirplex  conjugate  pairs.  Calling  the  real  and  imaginary 
parts  of  the  pairs  RE^ , AIM^,  RE,, , AIP^  equation  (15)  is  factored 
to  yield  two  cascaded  second  order  sections 


H(S)  = 


SWB 


SWB 


s - 2sre.  + re:  + aim: 


S - 2SRE-  + 


RE  ^ + 


aim: 


(16) 


For  each  0 of  a given  N,  the  program  calculates  roots  for  both 
sections  of  equation  (16)  and  labels  them  the  ith  and  the  ith  + 1 
section.  If  N,  the  number  of  second  order  sections  specified,  is 
even,  the  program  will  calculate  N pairs  of  RE  and  AIM  values 
or  2N  = n'  roots.  If  N is  odd,  the  last  value  of  0 is  0. 
Substituting  0=0  into  equation  (4)  and  factoring  yields  two 
identical  first  order  sections,  1/(S  + 1).  The  program  will  calculate 
N + 1 pairs  of  RE  and  AIM  values,  but  because  the  last  two  pairs 
are  the  same  due  to  the  identical  first  order  sections,  the  last  pair 
will  not  be  used. 

Because  both  second  order  sections  of  equation  (16)  are  of  the 
same  format,  we  will  deal  with  only  one  section,  the  ith  section 
and  let 


-2REi  = D 
RE^  + AIM 


1 

2 
i 


(17) 


The  design  of  an  n'th  order  band-pass  or  band-stop  filter  leads  to 
n'/2  second  order  sections.  Substituting  equations  (17)  into  one 
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section  of  equation  (16)  yields  the  transfer  function  for  the  ith 
section 


H (S)  = -S-^ 

S + SD.  + C. 

1 1 


The  extended  bilinear  Z transform,  equation  (2)  is  used  to  get  to 
the  digital  domain.  Employing  equation  (2)  on  equation  (18)  yields 
H^(Z)  for  the  ith  second  order  section. 


Hi(Z) 


2 2 2 
^WBZ  - ^WB 

~ , 2D~.  I 7 2d7 

z + — + V + z<2Ci  - ^ + ^ + ci> 


Putting  the  denominator  of  equation  (19)  in  monic  form  yields  the 
transfer  function  for  the  ith  second  order  stage  of  the  filter 


Ht(Z)  = 


KU(A0Z  + AXZ  + A2) 

z2  + BHz  + b2i 


This  equation  is  the  same  as  equation  (1)  with  the  exception  of  the 
su.o_.ipts.  For  all  four  filter  types  discussed  here,  the  scale 
factor,  , and  coefficients  and  B2  are  a function  of  the 

section  calculated,  while  the  coefficients  A^,  A^,  and  A2  are  the 
same  for  all  sections  calculated.  In  going  from  equation  (19)  to 
equation  (20)  we  have 


ST  i o)T 

Letting  Z = e = eJ  for  S = ju>  and  taking  the  magnitude  of 


H^jw)  we  have 


I «i ( jw)  I = Ku 


^AqCosCZiuT)  + A^cosfaT)  + k^)2  + (AgSin(2u>T)  + A^sinCwT))2 
/(cos(2u)T)  + Blicos(wT)  + B21>2  + (sin(2wT)  + Bu9in(u)T))2 


The  magnitude  function,  equation  (22)  is  the  same  for  all  the  filters 
discussed  in  this  report. 

C.  Butterworth  Band-Stop  Filter 

The  design  procedure  is  almost  exactly  the  same  as  that 
of  the  Butterworth  band-pass  filter,  except  that  the  transformation 
to  band-stop  is  the  reciprocal  of  equation  (14),  i.e. 
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and  we  find  11^(5)  to  be 


H.(S) 


2 2 
S + WDM 

S2  + SD.  + C. 

l l 


(24) 


After  employing  the  extended  bilinear  Z transform,  equation  (2), 
we  have 


A0  ■ *2  ' ^2  + “Dm2 

A = 2 WDM2  - 
1 T 


(25) 


and  B^ , B21,  ^ are  the  same  functions  of  and  as  in 

equation  (21) . These  coefficients  are  then  used  in  the  calculation 
of  equation  (22)  to  find  |H  (jw)|. 

D.  Chebychev  Band-Pass  Filter 

The  Chebychev  filter  ripples  with  equal  amplitude  in  the 
pass-band.  The  amount  of  ripple  is  specified  by  the  quantity  6 
(labeled  RIP  in  the  program).  The  poles  of  the  filter  are  found  on 
an  ellipse  described  by  two  Butterworth  circles  of  radii  A and  B 
"ith  A < B.  The  location  of  the  poles  on  the  ellipse  is  a function  of 
the  ripple  and  is  given  by  the  following  equation: 

B,  A = |-((*4"2  + 1 + e_1)1,/N  ± (Jz~2  + L + e-1)-1^)  (26) 


where 


.(1  - <s)2 


(27) 


and  N is  numerically  equal  to  the  order  of  the  low-pass  filter 
which  is  transformed  to  yield  the  band-pass  filter.  B is  given 
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for  the  plus  sign  and  A for  the  minus  sign.  The  Chebychev  ellipse 
then  has  major  axis  B and  minor  axis  A.  The  location  of  the  S 
plane  poles  on  the  ellipse  is  given  by 


Real  Part  = A cos0 
Imaginary  Part  = 3 sinO 


(28) 


The  0 ' s are  the  same  as  given  for  the  corresponding  order 
Butterworth  filter.  An  example  of  Chebychev  pole  locations  is 
illustrated  in  Figure  2.  For  A = y and  B **  1 in  a fourth  order 
filter,  0 = 22.5°  and  67.5°.  The  Chebychev  pole  locations  are 
determined  from  equations  (26),  (27)  and  (28). 

The  analog  second  order  Chebychev  low-pass  filter  is 


H(S) 


K2 

1 

2/N 

/ 2 

L/l  + e J 

S2  + KgS  K2 


(29) 


e is  calculated  from  equation  (27)  and  N is  equal  to  the  order  of 
the  low-pass  filter  which  is  transformed  to  yield  the  band-pass  filter. 


Kg  ind 

K2  are  calculated  by 

K_  = 2Acos0 

o 

(30) 

2 2 2 2 

K2  = A cos  0 + B sin  0 

• 

(31) 

The  substitution  of  the  low-pass 

to  band-pass  transformation. 

equation 

(14),  into  equation  (29) 

yields 

2 2 

S WB  K2 

1 

2/N 

A -t£2. 

S4  + S3KgWB  + S2(2WDM2  + K2WB2)  + SKgWDM2WB  + WDM4 


H(S) 


. (32) 


r 


f 
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After  finding  the  roots  of  equation  #32)  and  making  the  substitutions 
given  by  equations  (17)  we  find  the  ith  second  order  section 


SWK„ 


Vs)  ’ j 


S + SD.  + C, 
l i 


K3  - /q, 


fl  + 


1/N 


(33) 


Applying  the  extended  bilinear  Z transform  equation  (2)  yields  an 
equation  of  the  form  of  equation  (20)  where 


(34) 


B^  and  B2i  are  the  same  functions  of  and  given  by 

equations  (21).  These  coefficients  are  then  used  in  equation  (22)  to 
find  |Hi(ju))|. 


C.  Chebychev  Band-Stop  Filter 

Given  equation  (29)  for  H(S)  we  apply  the  low-pass  to 
band-stop  transformation  equation  (23)  to  obtain  the  4th  order 
transfer  function 


2 2 2 
(S*  + WDM  j n. 


\a7?. 


2/N 


H.  (S) 


j \'3>  a •>  2 2 2 2 

K2S  + S KgWB  + S (WB  + 2K2WDM  ) + SKgWDM  WB  + K2 


. (35) 


WDM 


N is  equal  to  the  order  of  the  low-pass  filter  which  is  transformed 
to  yield  the  band-stop  filter. 
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After  finding  the  roots  of  equation  (35)  and  making  the 
subi'.tit"*-1  on?  gi”°n  b”  ennarions  (17)  the  ith  second  order  section  is 


(S2  + wdm2)k. 


Hi(S)  = 


S + SDi  + Ci 


k3=  ^ 


1 + e 


1/N 


(36) 


Applying  the  extended  bilinear  Z transform  equation  (2)  yields  an 
equation  of  the  form  of  equation  (20)  where 


A0  ■ A2  ■ ^ + TOm2 


A = 2 WDM2  - 
1 T 


(37) 


Kli  G. 

l 


B.  . and  B„ . are  the  same  functions  of  C.  and  D given  by 
li  2i  i 1 


equations  (21).  These  coefficients  are  then  used  in  equation  (22) 


to  find  | (ju) | . 
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II . Using  the  Program 

The  first  data  card  read  into  the  program  contains  the  number  of 
second  order  sections  to  be  cascaded,  N,  and  the  type  of  filter 

desired,  KN.  N is  equal  to  1,  2 or  6,  which  corresponds  to 

the  order  of  the  low-pass  filter,  and  hence  corresponds  to  the  2nd, 
4th,...,  or  12th  order  band  pass  or  band  stop  filter  respectively. 

KN  is  the  type  of  filter  desired.  The  values  of  KN  specifies 


one  of  the  four  choices  given  by 


KN  Type 

1 Butterworth  Band-Pass 

2 Butterworth  Band-Stop 

3 Chebychev  Band -Pass 

4 Chebychev  Band-Stop 

The  format  on  the  N,  KN  card  is  212. 

The  second  data  card  read  in  is  the  sampling  interval  T in 

F10.6  format.  When  choosing  T,  1/T  should  be  approximately  equal 

to  ten  times  the  center  frequency  (WDM) . 

The  third  data  card  read  in  contains  the  values  of  the  upper 

and  lower  cutoff  frequencies,  and  u^,  in  2F10.4  format.  For 

the  Butterworth  filters,  the  cutoff  frequencies  are  the  -3db  cutoff 

frequencies.  For  the  Chebychev  filters,  the  magnitude  of  the 

2 x' 

response  is  1/(1  + e ) 2 = 1 - <$  at  the  cutoff  frequencies,  w is  in 
radians.  <S  is  the  ripple  factor. 

If  the  desired  filter  is  Chebychev,  i.e.,  KN  * 3 or  4,  the 
next  data  card  contains  the  ripple  (RIP)  factor  in  F5.3  format. 

If  the  desired  filter  is  Butterworth,  i.e.,  KN  = 1 or  2,  this 
card  is  omitted  from  the  data  deck. 
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The  final  data  card  is  the  starting  frequency  (FREQ1)  and  the 
frequency  increments  (DELT)  in  radians.  The  format  of  the  FREQ1, 
DELT  card  is  2F10.4.  Determine  DELT  by  the  following: 


DELT  = 


final  frequency  - starting  frequency 


This  is  necessary  because  there  are  1024  frequency  data  points 


calculated  in  the  program.  Choose  FREQ1  and  DELT  to  insure  that 


calculated  values  will  include  the  data  of  interest.  For  maximum 
efficiency  of  the  program,  DELT  should  be  a multiple  of  2 so  no 
decimal  to  binary  conversion  errors  are  incurred. 

The  digital  filter  coefficients  are  computed  and  printed  out  for 
each  second  order  section.  The  full  filter  magnitude  response,  as  well 
as  each  section  magnitude  response,  is  printed  for  each  of  the 


specified  frequency  increments.  When  there  is  only  one  second  order 
section,  the  section  magnitude  response  j.s  the  full  filter  magnitude 
response  and  is  only  printed  once. 

The  program  may  be  easily  modified  to  incorporate  a graphics 
display  of  the  magnitude  response.  There  is  a comment  card  in  the 
BPASS  program  indicating  where  the  graphics  subroutine  call  card 
should  be  inserted. 

The  program  is  written  with  input  obtained  via  device  4 and 
output  written  to  device  6.  These  numbers  should  be  assLgned  to  the 
appropriate  devices  prior  to  running  the  program. 

The  program  was  developed  on  a PDP-11/20  with  a DOS/BATCH 


operating  system.  Trial  runs  frequently  used  a TTY  terminal  as  well 
as  a card  reader  for  input  (device  4) ; and  a TTY  terminal  as  well 
as  a line  printer  for  output  (device  6).  Double  precision  arithmetic 
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is  employed.  To  decrease  required  memory  storage,  only  the  frequency 
interval  values  and  the  full  magnitude  response  are  saved.  The 
section  magnitude  responses  are  printed  out,  but  are  not  stored. 

The  program  will  produce  approximately  21  pages  of  output. 


Shown 

below  are  sample  deck 

set-ups. 

Data 

Card 

Format 

Fxample 

1 

212 

0504  (5  sectiors  Chebychev  banc-si  o) 

O 

c. 

F10.6 

0.002  (T  = 0.002) 

3 

2F10.4 

60  40  («u  = 60,  = 40  radians) 

4 

F5.3 

0.10  (Ripple  amplitude  = 0.10) 

S 

2F10.4 

0 0.1  (Start  at  u = 0.  Steps  of 

0.1  radian.  Will  finish 
just  past  a)  = 102  radianr . ) 

1 

212 

0401  (4  sections  Butterworth  band-n.'ss) 

■> 

F10.6 

0.002  (T  = 0.002) 

3 

2F10.4 

60  40  (u>u  = 60,  = 40  radians) 

4 

2F10.4 

0 0.1  (Start  at  to  *=  0.  Steps  of 

0.1  radian.  Will  finish 
just  past  u)  » 102  radians)  . 

The  following  pages  contain  annotated  examples  of  output  data. 


w H HI  H2  H3  H4 

0.000  0.10000E+01  0. 11726E+01  0.85284E+00  0.68611E+00  0.14575E+01 

0.1000  0 . 10000E+01  0 . 11726E+01  0.85284E+00  0.68611E+00  0.14575E+01 

0.2000  0. 10000E+01  0.11726E+01  0.85284E+00  0.68611E+00  0.14575E+01 

0.3000  0.99999E+00  0.11726E+01  0.85283E+00  0.68610E+00  0.14575E+01 

0.4090  0 . 10000E+01  0.11726E+01  0.85283E+00  0.68609E+00  0.14576E+01 

0.5000  0 . 10000E+01  0.11726E+01  0.85282E+00  0.68609E+00  0.14576E+01 


-171- 


WDU  is  the  prewarped  upper  frequency. 

WDL  is  the  prewarped  lower  frequency. 

WDM  is  the  prewarped  center  frequency. 

WB  is  the  bandwidth,  WDU  - WDL. 

T is  the  sampling  interval. 

The  Real  and  Imaginary  part  of  the  roots  of  the  filter  are  given  next. 

I is  the  ith  stage.  I varies  from  1 to  N. 

Aq,  , A2  are  the  Butterworth  band-stop  filter  numerator  coefficients, 
is  the  gain  factor. 

B^ , and  B.?  are  the  Butterworth  band-stop  filter  denominator  coefficients. 
W is  the  frequency 

H is  the  overall  magnitude  of  the  digital  transfer  function 
HI  is  the  magnitude  of  the  digital  transfer  function  (1st  stage). 

H2  is  the  magnitude  of  the  digital  transfer  function  (2nd  stage). 

H3  is  the  magnitude  of  the  digital  transfer  function  (3rd  stage). 

H4  is  the  magnitude  of  the  digital  transfer  function  (4th  stage). 


See  Figure  4. 
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This  is  an  example  of  the  output  for  an  8th  order  Chebychev 
band-pass  filter  (0403)  with  T = 0.002,  io^  = 60  radians,  and 

= 40  radians.  The  starting  frequency  is  0 radian,  the  frequency 
increment  is  0.1  radian,  and  the  ripple  is  0.1. 


WDU  = 60.07210  WDL  = 40.02135  WDM  = 49.02902  WB  = 20.05076 
T = 0.20000E-02 


A = 0.37642105 
A = 0.37642105 


B = 1.06850027 
B = 1.06850027 


K = 0.69553541 

K„  = 0.28810020 
o 


K = 0.28813942 
= 0.99524620 


THE  ROOTS  OF  THE  FILTER  ARE  GIVEN  BELOW 
REAL(l)  = -3.19528085  IMAGINARY (1)  = -44.97792290 

REAL(2)  = -3.77772492  IMAGINARY (2 ) = -53.17662810 

REAL (3)  = -1.15829660  IMAGINARY (3)  = -40.10115290 

REAL(4)  = -1.73001696  IMAGINARY (4)  = -59.89456990 


THE  COEFFICIENTS  OF  EACH  DIGITAL  FILTER  SECOND  ORDER  SECTION  ARE 
GIVEN  BELOW 


FOR 

I = 

1 

A° 

0 . 10000000E+01 

A1 

= 

0.00000000E+00 

-0. 10000000E+01 

K! 

= 

0. 10395603E-01 

= 

-0 . 19792607E+01 

B2 

= 

0.98732565E+00 

FOR 

I = 

2 

‘2 

a 

0 . 10000000E+01 

A, 

= 

0.00000000E+00 

= 

-0 . 10000000E+01 

K1 

= 

0.10375296E-01 

= 

-0. 19737935E+01 

B2 

= 

0 . 98504460E+00 

FOR 

I = 

3 

A„ 

s 

0. 10000000E+01 

at 

= 

0.00000000E+00 

fl2 

= 

-0 . 10000000E+01 

K1 

= 

0 . 19406846E-01 

= 

-0. 19889723E+01 

B2 

= 

0 . 99538493E+00 

FOR 

I = 

4 

A0 

= 

0 . 10000000E+01 

A1 

S 

O.OOOOOOOOE+OO 

*2 

= 

-0 . 1000P000E+01 

= 

0. 19346636E-01 

» 

= 

-0.19788675E+01 

B2 

= 

0 . 99312838E+00 

W 

H 

HI 

H2 

H3 

114 

0.0000 

0 . OOOOOE+OO 

0. 00000E+00 

0. 00000E+00 

0. 00000E+00 

0. 00000E+00 

0.1000 

0. 12493E-12 

0 . 51560E-03 

0. 36886E-03 

0.12106E-02 

0.54265E-03 

0.2000 

0. 19990E-11 

0. 10312F-02 

0.73773E-03 

0 . 24211E-02 

0.10853E-02 

0. 3000 

0. 10121E-10 

0. 15468E-02 

0.11066E-02 

0 . 36  318E-02 

0 . 16280E-02 

0.4000 

0. 31991E-10 

0. 20625E-02 

0 . 14755E-02 

0 . 48426E-02 

0.21707E-02 

n. 5000 

0.7811 6E-1 0 

0 . 25  783E-02 

0 . 1 8445E-02 

0.60536E-02 

0.27134E-02 

A 
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WDU  is  the  prewarped  upper  frequency. 
WDL  is  the  prewarped  lower  frequency. 
WDM  is  the  prewarped  center  frequency 
WB  is  the  bandwidth,  WDU  - WDL. 

T is  the  sampling  interval. 

B 


, A = !((£*  + 1 + e_1)1/N  ± (/e-2  + 1 + e-1)"1/N) 

Kg  = 2Acos(0) 

K2  = A2cos2(e)  + B2sin2(0) 

Ihe  Real  and  Imaginary  part  of  the  roots  of  the  filter  are  given  next. 

I is  the  ith  stage.  I varies  from  1 to  N. 

Aq,  A^ , A2  are  the  Chebychev  band-pass  filter  numerator  coefficients, 
is  the  gain  factor. 

B^,  and  are  the  Chebychev  band-pass  filter  denominator  coefficients. 
W is  the  frequency. 

H is  the  overall  magnitude  of  the  digital  transfer  function. 

HI  is  the  magnitude  of  the  digital  transfer  function  (1st  stage). 

H2  is  the  magnitude  of  the  digital  transfer  function  (2nd  stage). 

H3  is  the  magnitude  of  the  digital  transfer  function  (3rd  stage). 

H4  is  the  magnitude  of  the  digital  transfer  function  (4th  stage). 

See  Figure  5. 


Figure  2 
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MAGNITUDE  VS  FREQUENCY 
FOR 

DIGITAL  TRANSFER  FUNCTION 
8th  ORDER  BUTTERWORTH  BAND-PASS  FILTER 


N = 4 

Start  at  oj  = 0 radian 
T = 0.002 

Steps  of  0.1  radian 
oj  =60  radians 
= 40  radians 


oj  (radians) 


Figure  3 
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MAGNITUDE  VS  FREQUENCY 
FOR 

DIGITAL  TRANSFER  FUNCTION 
8th  ORDER  CHEBYCHEV  BAND-PASS  FILTER 


N = 4 

Start  at  w = 0 radian 
Ripple  =o=0. 100 
T = 0.002 

Steps  of  0.1  radian 
u.1  = 60  radians 

= 40  radians 


w (radians) 


Figure  5 
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Abstract 


The  use  of  a parameter  identification  procedure  to  detect 
faults  in  hardware  used  to  implement  a broad  class  of  linear 
algorithms  defined  as  digital  filters  is  presented.  Using  the 
filter  coefficient  estimates  produced  by  the  identifier  a method 
of  measuring  the  acceptability  of  the  filtering  algorithm  is 
suggested  and  a numerical  example  is  given. 
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INTRODUCTION 

Today's  inteqrated  circuit  technoloqy  has  provided  inexpensive  diqital 
hardware  that  can  be  used  for  the  direct  implementation  of  the  diqital 
signal  processing  algorithms  utilized  in  a variety  of  communication  and 
control  systems.  In  these  applications,  there  is  a need  for  new  methods 
of  failure  analysis.  While  it  is  desirable  to  know  when  a hardware  failure 
occurs,  and  to  know  where  the  failure  is  located,  it  is  also  important  to 
know  what  effect  the  failure  has  on  algorithm  performance.  In  many  cases 
a failure,  such  as  losing  the  least  significant  bit  in  a register,  will  not 
significantly  degrade  performance  and  there  is  no  need  to  consider  switching 
to  a redundant  implementation. 

Currently,  most  of  the  work  on  failure  analysis  is  described  under  the  heading 
of  fault  detection.  From  an  operational  point  of  view,  desiqners  are  concernec 
with  the  development  of  fault  tolerant  computer  systems.  However,  in  both 
cases  the  hardware  is  usually  considered  and  little  attention  is  given  to 
the  interaction  between  the  hardware  and  the  algorithm  to  be  implemented. 

To  provide  more  information  about  algorithm  performance,  the  problem  of 
mathematically  describing  failure  detection  and  diagnosis  has  recently  been 
given  attention.  Mehra  and  Peschon  [1]  suggested  implementing  a Kalman 
filter  in  parallel  with  an  algorithm  implementation  and  usinq  the  statistics 
of  the  innovations  process  for  detecting  system  failures.  Davis  [2]  shows 
a method  for  usinq  a Kalman  filter  to  estimate  the  time  when  a failure 
occurs.  He  then  recommends  readjustment  of  the  filter  parameters  to  ODtair 
new  state  estimates  after  the  occurence  of  the  fault.  Krischer  [3],  in  an 
applications  oriented  approach,  applied  a parameter  identification  approach 
t.o  estimate  the  state  of  a biological  system  as  the  system  parameters  slow  . 


vary. 


-182- 


In  this  paper,  a parameter  identification  procedure  similar  to  that  presc.'.tc- 
by  Mendel  and  Fu  T4]  is  used  to  detect  faults  in  the  implementation  of  a 
broad  class  of  algorithms  defined  as  digital  filters.  Given  a design,  the 
filter  coefficients  are  used  as  initial  conditions  for  the  identifier  operating 
in  parallel  with  the  filter.  The  identifier  output  consists  of  an  estimate 
of  the  filter  coefficients  and  an  error  signal.  When  a fault  occurs,  the 
error  signal  and  the  coefficient  estimates  will  change  rapidly  during  the 
next  few  sampling  times.  If  the  fault  is  very  serious,  such  as  complete 
failure  of  the  multiplier,  the  algorithm  implementation  no  longer  exists  anc 
a total  failure  has  occured.  If,  however,  the  failure  is  not  total  so  that 
the  effect  is  to  modify  the  filter  coefficients,  the  coefficient  estimates 
will  converge  to  the  new  values.  At  this  point,  the  effect  of  the  coefficient 
changes  must  be  evaluated  to  determine  if  the  algorithm  characteristics  are 
still  satisfactory.  In  this  paper  this  is  done  by  establishing  bounds  on 
the  steady-state  magnitude  and  phase  functions.  If  the  change  in  the 
coefficients  allows  these  bounds  to  be  satisfied,  the  filter  operation  is 
considered  to  be  acceptable. 

The  procedures  described  have  several  advantages.  First,  remote  sensing  on 
an  alaorithm  implementation  can  be  accomplished.  Secondly,  the  identifier 
can  be  used  with  any  linear  system  that  can  be  modeled  by  a difference 
equation.  Thirdly,  a better  assessment  as  to  the  need  for  switching  in  an 
alternate  algorithm  implementation  is  available. 
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THfc  DIGITAL  FILTER  MODEL 


An  nth  order  time  invariant  linear  digital  filter  is  represented  by  the 
expression 


y[nT'i=£  a.  x [(n-K)T]  - £ t>,y[(n-j)T] 

k=G  k M 3 


(1) 


Where  the  coefficient  sequences  {a^}  and  { b ^ } are  chosen  to  uchie ' rn. 
desired  filter  characteristics.  Usinq  a vector  representation  (1)  can  be 
rewritten  as 


y[nT  1 = a^fnT]  - bV  [(n-l)T' 


wnere 


a [aQ  a^  . • am]» 
-t  = [bl  b2  • * bj. 


(2) 

(3) 

(4) 


xfnT]  = 


[x[nT] 

x[(n-l)T] 

xf (n-m)T] 


(5) 


and 


y [ (n-'l  )T]  = 


y[(n-l)T 

y[(n-2)T 

ly[0T] 


(6) 


Letting 


(7) 


zfnT j = y[nT] 


(P) 
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and 


u[nT]  = 


x[nT] 

XL(n-l)T] 


(2)  can  be  rewritten  as 


(9) 


z[nTl  = ^[nTjc  = <u[nT],c>  (1C; 

Equation  (10)  now  represents  the  model  for  the  digital  filter  that  is  use^ 
in  me  identification  algorithm. 


THE  IDENTIFICATION  ALGORITHM 

The  identification  procedure  is  a sequential  gradient  descent  algorithm 
utilizing  a quadratic  cost  function.  The  model  is  depicted  in  the  block 
diagram  of  Fig.  1.  Here,  at  time  t=nT  the  actual  digital  filter  output  is 
y[nT]  and  the  current  input  is  jj[nT].  From  (10)  the  filter  model  is  now 
expressed  as 

z[nT]=kt[nT]c[nT]  (11) 

Where  £[nT]  is  the  approximation  of  the  parameter  vector  c at  time  t=nT. 
Defining  the  error  as 

e[nT]-y[nT]-z[nT]  (12) 

the  quadradic  error  function  is 

J[c[nT]>4[y[nT]-z[nT]]2  (13) 

To  minimize  J[c[nT]]  a multidimensional  form  of  Newton's  method  is  used 
giving  the  recursive  expression 

c.[(n+l  )T]=c[nT]  <nT]  grad  [J[c[nT]]].  (14) 

Where  R[nT‘  is  an  nxn  matrix  whose  coefficients  are  yet  to  be  specified. 
From  (13)  the  expression 
graa  JFc[nT]]=-eLnT]u[nT] 

is  obtained.  Substitution  of  (15)  into  (14)  now  yields  the  recursive 


(15) 


-185- 


expression 

c[(n+l )T]=c[nT]+R[nT]e[nT]u£nT] 


Mendel  and  Fu  [7],  show  that  for  the  choice 

I 


R[nT]=-r 

£ LnT]u[nT] 

convergence  in  the  mean  for  clnT]  as  n--  °°  is  obtained. 


(16) 


(17) 


The  addition  of  the  measurement  noise  v[nT]  and  n[nT]  as  shown  in  Fig.  1 
causes  the  parameter  estimates  using  (16)  and  (17)  to  become  biased.  To 
avoid  this  a new  error  function 

e [nT]  = *?[y[nT]  - [z[nT]+v[nT]]]2  (18) 

is  used  to  give 

c[(n+1)T]  = [I+R[nT]$[nT]]c[nT]+R[nT]e[nT][£[nT]+n[nT]]  (T9) 

where  $[nT]  is  the  nth  estimate  of  the  measurement  noise  covariance. 


FILTER  PERFORMANCE 

Oiven  a digital  filter  described  by  (1),  the  implementation  of  the  filter 
is  usually  done  to  minimize  the  effects  of  coefficient  and  rounding  errors. 

In  Dractice  this  results  in  most  filters  being  implemented  as  a cascade 
or  parallel  connection  of  first  and  second  order  sections  [5]. 

Oiven  a cascade  or  Darallel  connection,  several  options  for  utilizing 
the  identifier  are  possible.  First,  the  identifier  can  treat  the  complete 
structure  as  a single  filter.  However,  for  hiah  order  filters,  the  deter- 
mination of  the  coefficient  sensitivities  with  respect  to  any  error  criterion 
is  oenerally  very  difficult.  Secondly,  the  identifier  can  operate  on 
each  first  or  second  order  section.  This  allows  a performance  evaluation 
at.  the  section  level  which  can  be  interpreted  in  terms  of  overall  performance, 
because  magnitude  and  phase  are  often  used  to  specify  a digital  filter 
desinn,  a change  in  the  filter  coefficients,  due  to  a fault,  will  " e 
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evaluated  co  see  if  the  specifications  are  still  satisfied,  this  follows 
the  coefficient  desiqn  procedure  described  by  Brubaker  [6], 

For  a second  order  filter  with  the  transfer  function 


H(Z)  - 


Val7‘1+a2Z"2 

ri 

l+b,Z  '+b2Z 


(20) 


the  magnitude  function  is  described  by 


where 

A(u>)  = jaQ7+a.J'"+a?  >(2a0a1+2a1a2)  cosu>T  + 2aQa2  cos  2wt|  ^ (22) 

and 

3(u)  = |l  + b,?  + b22  + 2b1  (1  + b2)  cosuT  + 2b?  cos  2UT j 2 (23) 

Using  (21)  a differential  magnitude  approximation  can  be  written  as 


-ft 


.j.  »MMi 

aao 


Aa, 


alH(jui) 

aa1 


Aa, 


a |H(,iu3) 

3an 


Aa, 


" | H ( ioj)  1 
3b, 


*»,  ♦ “>2 


(24) 


'he  partial  derivatives  in  (24)  are  given  in  [6]  and  by  evaluating  these 
.er  /atives  over  a frequency  ranqe  of  interest  a region  in  parameter  space 
a',  oe  established  where  filter  performance  is  satisfactory  for  a given 
A rt[jai!j.  A similar  strategy  can  be  Implemented  for  the  phase  function. 


expression 

A, 


jc[(n+l )T]=c[nT]+R[nT]e[nT]uJnT] 

Mendel  and  Fu  [7],  show  that  for  the  choice 

I 


R[nT]=— r 

£ [nT]u[nT] 

convergence  in  the  mean  for  c[nT]  as  n->-  ® is  obtained. 


06) 


(17) 


The  addition  of  the  measurement  noise  v[nT]  and  n[nT]  as  shown  in  Fig.  1 
causes  the  parameter  estimates  using  (16)  and  (17)  to  become  biased.  To 
avoid  this  a new  error  function 

e [nT]  * *s[y[nT]  - [z[nT]+v[nT]]]2  (18) 

is  used  to  give 

c[(n+l)T]  = [I+R[nT]$[nT]]c[nT]+R[nT]e[nT]i;u[nT]+n[nT]]  (19) 

where  $[nT]  is  the  nth  estimate  of  the  measurement  noise  covariance. 


FILTER  PERFORMANCE 

Given  a digital  filter  described  by  (1),  the  implementation  of  the  filter 
is  usually  done  to  minimize  the  effects  of  coefficient  and  rounding  errors. 

In  oractice  this  results  in  most  filters  being  implemented  as  a cascade 
or  parallel  connection  of  first  and  second  order  sections  [5]. 

Given  a cascade  or  oarallel  connection,  several  options  for  utilizing 
the  identifier  are  possible.  First,  the  identifier  can  treat  the  complete 
structure  as  a sinqle  filter.  However,  for  hiah  order  filters,  the  deter- 
mination of  the  coefficient  sensitivities  with  respect  to  any  error  criterion 
is  oenerally  very  difficult.  Secondly,  the  identifier  can  operate  on 
each  first  or  second  order  section.  This  allows  a performance  evaluation 
at.  the  section  level  which  can  be  interpreted  in  terms  of  overall  performance, 
because  magnitude  and  phase  are  often  used  to  specify  a digital  filter 
design,  a change  in  the  filter  coefficients,  due  to  a fault,  will  be 
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EXAMPLE 

To  illustrate  using  identification  for  fault  detection  consider  the  first 
order  digital  filter 


H(Z)  = 


V 


-1 


0.8Z 


-1 


1-b.jZ 


1 -0.85Z* 


(15) 


For  A i H C.ito]  1 = 0.1,  the  acceptable  region  for  the  filter  coefficients  is 
shown  in  Fig.  2.  The  identifier  response  to  the  filter  is  shown  in  Fig.  3 
ds  b^  chanaes  from  0.85  to  0.4.  Initially  the  identifier  tracks  the 
correct  coefficients.  When  bj  changes  the  error  signal  changes  rapidly 
and  converges  to  zero.  The  coefficient  estimates  converge  to  the  new 
values  and  the  b^  coefficient  of  0.4  does  not  allow  satisfactory 
oerformance  through  use  of  the  region  shown  in  Fig.  2.  A redundant  filter 
would  then  be  set  into  operation. 
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Figuro  1.  Block  Diagram  for  the  Gradient  Identifier 


Ar’i)  l i Hide 
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